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1.  INTRODUCTION 


The  aim  of  this  report  is  to  present  new  elements  in  a  rigorous  least- 
squares  adjustment  of  geodetic  quantities.  Although  the  initial  emphasis  was  on 
the  satellite-to  satellite  tracking  (SST),  the  adjustment  aspects  developed 
during  the  course  of  the  study  are  broad  enough  to  be  applicable  to  a  variety  of 
problems  in  physical  sciences.  The  body  of  the  report  is  thus  reserved  for  the 
nonlinear  least-squares  method  with  or  without  constraints,  whereas  the  SST 
model  and  its  adjustment  are  described  in  Appendix  A.  The  constraints  are 
considered  to  be  in  general  nonlinear;  linear  constraints  Joined  to  a  linear 
model  are  presented  merely  as  a  special  case.  A  relatively  new  category  of 
constraints,  termed  inequality  constraints,  is  treated  in  Appendix  B. 

The  SST  model,  when  considered  without  any  approximations  or  deformations, 
such  as  the  smoothing  and  differentiation  of  the  original  observations,  appears 
to  have  the  form  of  the  general  adjustment  model.  The  latter  usually  features 
the  observables  and  the  parameters  interwoven  in  nonlinear  relationships. 
However,  the  observables  in  the  SST  model  have  special  characteristics  allowing 
for  a  substantial  simplification  In  partlculaj  .  they  are  present  in  linear 
combinations  involving  merely  two  observables  per  equation.  A  simple  linear 
transformation  of  the  observables,  which  affects  the  original  observations,  the 
adjusted  observations,  and  the  residuals  in  the  same  fashion,  changes  the 
general  model  into  the  parametric  adjustment  model.  The  latter  is  treated  here 
In  a  nonlinear  form.  One  can  thus  avoid  the  approximations  caused  by  the 
truncation  of  second-  and  higher  order  terms  in  the  model's  Taylor-series 
expansion,  as  is  done  In  the  standard  linearized  approach.  The  transformation 
of  observables  leading  to  the  parametric  model  can  be  useful  not  only  in  view  of 
the  rigorous  least-squares  adjustment  of  the  SST,  but  in  the  treatment  of  other 
kinds  of  geodetic  data  as  well. 

The  above  transformation  has  allowed  the  analysis  to  shift  focus  from  the 
general  adjustment  model,  linear  or  nonlinear,  to  the  nonlinear  parametric 
adjustment  model.  The  linear  version  of  the  latter  follows  as  a  simple  special 
case;  it  has  been  thoroughly  analyzed  and  extensively  used  over  the  past  several 
decades.  The  least  squares  adjustment  of  the  nonlinear  parametric  model  is  the 
subject  of  Cha|)ter  2,  summarizing  the  recent  development  contained  in  the  AFGL 
report  [Blaha,  1989J .  This  development  is  based  on  analogies  between 


adjustments  and  geometry  which  have  led  to  the  conception  of  an  isomorphic 
geometricaJ  setup. 

In  some  adjustment  problems  the  parametric  model  is  subject  to  constraints, 

l.e.,  conditions  to  be  fulfilled  by  its  parameters.  As  an  illustration  related 

to  the  SST,  one  notices  that  although  equation  (2)  in  Appendix  A  does  not 

contain  any  constraints  among  the  parameters,  ('ertain  types  of  constraints  may 

be  needed  in  conjunction  with  the  geopotontial  representation  used  by  the  model, 

and  with  other  factors.  In  an  early  treatment  of  the  SST.  Schwarz  (1970] 

considers  the  gravity  field  split  into  the  reference  field  described  by  an  (N.N) 

spherical -harmonic  expansion,  and  the  residual  field  described  by  density-layer 

parameters.  Because  of  the  initial  stipulation  that  the  density  layer  solution 

2 

should  not  affect  any  of  the  (N*-!)  coeff i<rients  in  the  underlying  spherlcal- 

2 

harmonic  model,  he  concludes  that  (N<^1)  constraints  should  be  included,  at 
least  in  theory.  Clearly,  this  notion  is  not  tied  only  to  the  density  layer 
parameters  since  similar  reasoning  applies  also  with  regard  to  other  localized 
representations  of  the  residual  gravity  field. 

In  analogy  to  an  adjustment  model  itself,  a  set  of  constraints  can  also  be 
linear  or  nonlinear.  A  constrained  least -squares  adjustment:,  where  the  emphasis 
is  on  nonlinear  constraints,  is  developed  in  Chapter  3  herein.  This  endeavor  is 
again  based  on  the  isomorphism  between  adjustments  and  geometry.  The  geometric 
derivations  leading  to  the  final  adjustment  formulas  are  carried  out  with  the 
aid  of  tensor  structure  and  notation.  It  should  be  mentioned  that  all  the 
constraints  considered  so  far,  whether  linear  or  nonlinear,  are  the  familiar 
equality  constraints. 

However,  recent  years  have  witnessed  an  increasing  interest  in  applications 
of  linear  inequality  constraints,  which  can  be  used  with  advantage  in  the 
problems  where  the  smoothing  effect  of  the  standard  least-squares  method  is 
undesirable.  The  concept  of  inequality  constraints  is  instrumental  in  reducing 
the  solution  space  to  a  band,  as  is  descrilted  by  Frltsch  [1987]  in  conjunction 
with  a  linear  parametric  adjustment  model.  This,  in  turn,  allows  one  to 
minimize  the  maximum  error  and  thus  to  accommodate  the  worst  case.  Appendix  B 
herein  presents  the  least-squares  algorithm  applicable  to  linear  models  with 
linear  inequality  constraints,  which  is  developed  using  again  an  isomorphic 
geometrical  setup  with  tensor  structure  and  notation. 


2.  NONLINEAR  PARAMETRIC  LEAST-SQUARES  ADJUSTMENT 

2 . 1  Mathematical  Background 

The  parametric  adjustment  model  expresses  each  of  the  observables  in  terms 
of  parameters,  where  the  structure  linking  the  two  groups  of  variables  is,  in 
general,  nonlinear.  The  number  of  observables  is  denoted  by  n  and  the  number  of 
parameters  by  u,  where  n  must  be  greater  than  u  for  an  adjustment  to  take  place. 
The  adjustment  model  is  written  as 

I.®  -  F(X®)  , 

where  L**  and  X**  are  the  sets  (column-vectors)  of  adjusted  observations  and 
adjusted  parameters,  respectively.  This  chapter  describes  the  resolution  of  a 
nonlinear  model  through  an  isomorphic  geometrical  setup  with  tensor  structure 
and  notation.  Such  efforts  date  back  to  (Blaha,  1984],  which  treats  a  linear  or 
linearized  adjustment  model.  Later  papers  and  reports,  such  as  [Blaha,  1987], 
contain  an  initial  analysis  of  a  nonlinear  model.  The  most  recent  development 
in  this  area  is  described  in  [Blaha,  1989]. 

In  a  standard  adjustment  approach,  a  nonlinear  adjustment  model  is  subject 
to  the  Taylor  series  expansion  based  on  an  initial  set  of  parametric  values,  X°. 
The  terms  in  the  second  and  higher  powers  of  the  parametric  corrections  are 
neglected,  resulting  in  the  familiar  (linearized)  observation  equations.  In 
matrix  notation,  the  latter  are  expressed  by 

V  =  A  X  +  I,  , 

where  A  is  the  design  matrix,  X=X  -X  is  the  column- vector  of  parametric 
corrections,  V-L^-L^  is  the  column  vector  of  residuals,  and  L  =  L°-L^  is  the 
column-vector  of  constant  terms,  with  L°=F(X°)  representing  observables 
consistent  with  the  initial  set  of  parameters,  and  containing  the  actual 
observations.  The  linearized  model  is  subjected  to  the  least-squares  criterion 

T 

V  P  V  =  minimum  . 


where  P  is  the  weight  matrix  of  observations.  This  criterion  leads  to  the 
formation  of  the  familiar  normal  equations. 


If  the  original  adjustment  model  is  nonlinear,  the  resolution  of  tht 

linearized  model  does  not  yield  the  final  answers.  The  process  is  usually 

repeated  with  new.  updated  parameters  and  the  corresponding  changes  in  A  and  L. 

However,  the  variance-covariance  matrix  of  observations,  I,  as  well  as  the 

weight  matrix  P.  adopted  as  P=E  are  constant.  Thus,  the  matrix  of  normal 
T 

equations.  N=A  PA.  changes  only  due  to  A.  and  the  column -vector  representing 

T 

the  right-hand  side  of  normal  equations,  U-  A  PL.  changes  only  due  to  A  and  L. 
The  computation  of  the  updated  parametric  values  through  a  new  X  requires  the 
formation  and  the  inversion  of  a  new  N  in  each  iteration,  or  a  mathematically 
equivalent  procedure.  When  X  becomes  sufficiently  close  to  zero  the  Iterative 
process  is  terminated.  As  its  by-product,  the  latest  matrix  N  ^  is  adopted  as 
the  variance-covariance  matrix  of  adjusted  parameters. 

The  functional  relationship  between  the  observables  and  the  parameters 
lends  itself  to  a  geometrical  interpretation  and  treatment  Involving  spaces  and 
surfaces  generalized  to  higher  dimensions.  In  particular,  the  formulation 

X  =  x  (u  )  ,  r  1,2 . n  ,  a  =  1,2 . u  , 

representing  the  parametric  adjustment  model,  can  be  linked  to  the  Gauss  form  of 

I* 

a  surface  in  relation  to  the  surrounding  space,  where  x  are  the  space 
coordinates  and  u®  are  the  surface  coordinates.  The  Gauss  form  of  a  two- 
dimensional  surface  (u=2)  embedded  in  a  three-dimensional  flat  space  (n=3)  is 
described,  together  with  two  other  forms,  in  Chapter  6  of  [Hotine,  1969],  In 
[Blaha,  1984],  both  the  n  dimensional  "observational"  space  and  the  u 
dimensional  "model"  surface  were  considered  flat.  The  latter  was  thus  In 
reality  a  hyperplane.  Although  the  model  surface  is  now  intrinsically  a  curved 
space,  the  surrounding  space  is  again  flat,  and,  as  is  shown  below,  its 
coordinate  system  Is  characterized  by  a  constant  metric  tensor. 

f 

In  denoting  the  n  observables  by  x  ,  r  =  1,2 . ii,  anti  the  u  unknown 

parameters  by  u®,  a =1,2 . ii,  we  can  represent  a  nonlinear  parametric 

adjustment  model  by 

r  r,  a,  r  .r  ,  a  .  a  .  0 

x  =  X  (u  )  =  X  *-  A  Au  •  (1/2|  0  .  Au  Au 
o  a  tip 

+  (1/6)  Au®  Au^  Au^  ♦  .  .  .  .  (la) 

otpr 

.  Cl  a  at  , ,  i.  V 

Au  =  u  -  u  ,  (lb) 
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(X  T  V  (X 

where  represents  an  Initial  set  of  paraaeters  and  x^=x  (u^)  represents 
the  observables  consistent  with  this  set.  The  lower-case  Roaan  indices  range 
froa  1  to  n,  and  the  lower-case  Greek  Indices  range  from  1  to  u.  Tensor 
syabolism  implies  the  summation  convention  over  the  dummy  (repeating)  Indices. 

In  the  geometrical  context,  the  first  equality  in  (la)  represents  the  Gauss 
fora  of  a  u  dimensional  surface  eabedded  in  an  n  dimensional  space.  The  surface 
is  endowed  with  the  coordinate  system  {u*},  a=l,2 . u  and  is  referred  to  as 

p 

the  model  surface,  and  the  space  is  endowed  with  the  coordinate  system  (x  ), 
r=1.2 . n.  and  is  referred  to  as  the  observational  space.  The  second  equality 

p 

in  (la)  is  the  Taylor  series  expansion  of  x  from  the  "initial"  point  P  lying  in 

the  model  surface,  whose  model -surface  coordinates  are  u*  and  whose 

o 

p 

observational -space  coordinates  are  x^.  The  notation  identifying  the  partial 

derivatives  at  P,  such  as  cl  x^/ a  u*  =a'' ,  d^x^/dxi'^du^so’^  etc.,  is  self- 

a  0/9 

evident.  The  actual  observations  can  be  thought  of  as  describing  the  point  Q  in 
the  observational  space,  which,  due  to  measuring  errors,  does  not  lie  in  the 
known  model  surface.  The  task  at  hand  consists  in  determining,  from  the 
observed  point  Q,  a  model-surface  point  satisfying  the  least-squares  criterion. 

In  the  adjustment  context,  the  variance-covariance  and  the  weight  matrices 
of  observations  depend  on  the  quality  of  measurements.  They  are  independent  of 
the  adjustment  model,  of  the  Initial  set  of  parameters,  of  the  outcome  of 
observations,  etc.  Thus,  for  a  given  observational  mode  they  are  constant.  In 
the  " tradi tiona ' "  identification  of  (Blaha,  1184],  variance-covariance  matrices 
correspond  to  associated  metric  tensors,  and  weight  matrices  correspond  to 
metric  tensors.  Accordingly,  we  represent  the  variance-covariance  matrix  of 

p  g 

observations  by  the  observat ional -space  associated  metric  tensor  g  ,  and  the 
weight  matrix  of  observations  by  the  observational -space  metric  tensor  g^^,  and 
state  that  both  tensors  are  independent  of  the  form  of  the  model  surface,  of  the 
initial  point  P.  of  the  observed  point  Q.  etc.,  leading  to  the  simplification 

g  ^  constant  .  (2) 

sr 

rs 

One  could  also  attribute  the  tensors  g  and  g^^  to  the  point  Q  and  state  that 
the  geometrical  setup  must  account  for  t)  located  anywhere  in  the  observational 
space.  In  turn,  (2)  implies  that  the  observational  space  must  be  flat. 


If  the  set  Sx  denotes  the  coordinate  dillerences  between  the  observed 
point  Q  and  the  desired  •odel -surface  point  denoted  P,  It  corresponds  to  the 
negative  residuals,  and  the  least  squares  criterion  corresijnnds  to 


i-2  .8  ,-r 

5s  ==  5x  g  3x 
sr 


miniBum  . 


The  quadratic  torn  (3)  represents  the  square  of  the  distance  between  Q  and  P. 
Therefore,  the  desired  "least-squares"  point  P  «ust  be  the  foot -point  of  the 
straight  line  dropped  orthogonally  from  Q  onto  the  model  surface.  We  note  that 
if  any  other  adjustment  condition  were  used  in  lieu  of  the  least-squares 
crite'lon,  the  minimum  distance  property  (3)  would  not  exist  and  the  geometric- 
tensorial  treatment  of  the  adjustment  theory  would  probably  be  much  more  complex 
if  not  impossible. 


2 . 2  Sumniai  y  of  the  Geometrical  Developmeitt 

A  convenient  approach  for  resolving  nonlinear  least  squares  problems 

consists  in  using  an  isomorphic  geometrical  setup  with  tensor  structure  and 

notation.  Such  a  link  is  highlighted  by  the  consideration  that  the  least 

squares  criterion  gives  rise  to  a  minimum  distance  property.  Among  the  basic 

correspondences,  the  number  of  observations,  n,  and  the  number  of  parameters,  u. 

define  the  dimensionality  of  the  observational  space  and  of  the  model  surface. 

respectively.  Since  the  constant  variance  covariance  matrix  of  observations,  L. 

rs 

corresponds  to  the  associated  metric  tensor  g  '  ,  and  the  weight  matrix  of 

observations,  adopted  as  L  * .  corresponds  to  the  metric  tensor  g  .  the 

sr 

f' 

observational  space  is  endowed  with  a  coordinate  system  {x  )  such  that 


g  -  constant  , 
sr 


g  =  constant 


b  r 

The  set  L  of  actual  observations  corresponds  to  the  set  x^  of  observational 

space  coordinates  describing  the  point  Q  All  possible  sets  of  adjusted 

observations  (subject  to  no  criterion)  correspond  to  the  Gauss  form  of  the  model 

surface  endowed  with  a  coordinate  system  (u®) 


r  r,  a, 

X  =  X  (u  )  , 


a  =  1.2, 


The  final  set  of  adjusted  parameters  X  .  corresponds  to  a  particular  set 

CK 

11  of  model  surface  coordinates  describing  the  h:.st  squares  point  P  The  set 


of  initial  parameters,  X^,  corresponds  to  the  set  u*  of  model -surface 

o 

coordinates  describing  the  initial  point  P.  The  final  set  of  parametric 

corrections,  X,  then  corresponds  to  these;  quantities  are  assumed 

o 

to  be  small  (termed  first  order)  The  final  set  of  adjusted  observations, 

h^=F(X^),  corresponds  to  a  particular  set  of  observational -space 

coordinates  describing  the  least -squares  point  P.  The  initial  point  P  Is 

r  r  ot 

described  by  these  coordinates  as  x  (u^),  reflecting  its  counterpart 

I-°"F(X°).  fhe  set  of  negative  constant  terms,  -  L  ^  -  F(X°)  ,  corresponds  to  the 

r  r  r 

contravariant  vector  <5x  -x,  -x  ,  while  the  set  of  negative  residuals, 

b  ^  °  - 

V-L  -L  .  corresponds  to  the  contravariant  vector  Ax  ^x^-x  .  The  initial 

design  matrix,  A,  which  in  standard  observation  equations,  V-AX-L,  relates  the 

parametric  corrections  to  the  residuals,  corresponds  to  the  design  tensor 
r  r  a 

A  =ax  .^du  evaluated  at  P.  On  the  other  hand,  the  standard  adjustment 

r  r 

approach  does  not  have  equivalents  of  Q  „  and  ®  „  ,  which  form  three-  and 

a^r 

four  dimensional  arrays,  respectively,  and  contain  the  second-  and  the  third- 

r  (X 

order  partial  derivatives  of  x  with  respect  to  u  ,  evaluated  at  P. 

The  geometrical  approach  is  based  on  a  direct  exploitation  of  the  relation 

a;  g^^  Ax*'  -  0  .  (4) 

.  s 

where  A„  represents  the  design  tensor  evaluated  at  the  least-squares  point 

p 

P,  and  equation  (4)  itself  represents  the  orthogonality  condition  at  P.  The 
outcome  of  the  geometrical  development  is  considered  in  two  methods,  called 
geometrical  and  extended  geometrical.  It  is  contrasted  to  the  standard  method 
treating  nonlinear  models  in  a  linearized  form.  The  algorithms  associated  with 
all  three  methods  are  presented  below  in  the  form  of  the  first  iteration,  and  in 
the  form  of  the  second  and  subsequent  iterations. 

In  tensor  notation,  the  initial  matrix  of  normal  equations  corresponds  to 
the  model -surface  metric  tensor  a^^  at  the  initial  point  P.  and  the  initial 
right-hand  side  of  normal  equations  corresponds  to  the  model-surface  covariant 

3 

vector  A-Ax  at  P,  where 
^  s 

a.  =  A®  g  A^  ,  Ax  -  g  Ax*^  . 

/3a  /9  sr  a  s  sr 


The  parametric  corrections  stemming  from  the  first  Iteration  are  symbolized  by 
(Au**),  and  they  give  rise  to  an  updated  point  (P).  The  latter  is  described  by 
the  mode  1 -surface  coordinates  ( u®  )  = u® ♦  ( Au® )  .  The  quantities  belonging  to  (P) 


will  likewise  be  written  in  parentheses.  The  parametric  corrections  obtained  in 
the  second  Iteration  will  be  denoted  A(Au®),  and  they  will  give  rise  to  a  new 
updated  point  determined  via  (u^" )  + A{Au®) .  The  notation  used  in  conjunction  with 
the  second  iteration  will  be  retained  also  for  any  further  iterations. 

<x  r 

standard  method.  Under  the  assumption  that  both  sets  Au  and  Ax  contain 
small  quantities  (first-order),  the  first  Iteration  in  the  standard  method  reads 

V 

representing  the  Initial  normal  equations.  The  second  and  further  iterations 
follow  the  same  principle: 

(a^^)  A(Au®)  =  (A‘^)  (6x^)  .  (5b) 

Geometrical  method.  Under  the  same  assumption  as  above  (both  sets  Au®  and 
Ax  contain  small  quantities),  the  first  iteration  utilizes  the  same  formula  as 
its  standard  counterpart,  namely 

(Au®)  =  A®  5k^  (6a) 

However,  the  second  and  further  iterations  proceed  according  to 

[(a^^)  -  (Ax^)  A(Au®)  =  (A^)  (Ax^)  .  (6b) 

representing  the  modified  normal  equations.  The  triply-indexed  quantity  0-  , 

po 

formed  by  second-order  partial  derivatives  of  the  observables  with  respect  to 
the  parameters,  is  evaluated  only  at  the  initial  point  P 

Extended  geometrical  method.  Although  the  assumption  regarding  Au®  is 
unchanged,  this  method  is  tailored  for  Ax  containing  relatively  large 
quantities,  for  which  t,*’e  first  Iteration  reads 

"la*  '  "e 

Compared  to  the  geometrical  method,  the  current  algorithm  is  seen  to  utilize 
second-order  partial  derivatives  and  to  giv'-  rise  to  the  modified  normal 
equatlc  .  already  in  Its  first  iteration.  The  formula  for  the  second  and 
furtr  I  teratlons  is  given  as 

'  S.'  *lar  ■  '''e> 
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representing  the  nodi  fled  noraal  equations  at  updated  stages.  Here  use  Is  nade 
s 


of  containing  third-order  partial  derivatives, 

evaluated  only  at  P,  siailar  in  this  respect  to  0®  . 

p« 


This  quantity  is 
We  note  that  the 


quantity  inside  the  brackets  of  (7b)  could  be  replaced  by  (a.  )-(3x  )(n-  ), 

pa  8  pa 

where  (0®  )  would  represent  the  second-order  partial  derivatives  evaluated  at 
pa 

an  updated  point. 


The  standard  adjustment  algorithm,  represented  by  the  relations  (5a, b) 
above,  results  in  the  projection  of  the  point  Q  onto  the  model  plane  passing 
through  the  initial  point  P,  followed  by  the  projection  of  Q  onto  a  new  model 
plane  passing  through  an  updated  point  (P).  etc.  The  orthogonality  condition 
(4)  is  then  fulfilled  essentially  as  a  by-product  of  these  projections.  By 
contrast,  the  geometrical  approach  actively  seeks  to  fulfill  it  at  every  step. 

A  one  -step  solution  producing  the  least-squares  point  P  directly  is  hindered 
only  by  the  necessity  to  truncate  certain  terms,  but  not  to  the  extent  of  making 
the  entire  model  linear  (see  the  above  equations  6b  and  7a, b).  The  matrix  of 
modified  normal  equations  generated  in  the  process  is  positive-definite,  similar 
in  this  respect  to  the  matrix  of  normal  equations  in  the  standard  method. 


Encouraging  results  have  been  obtained  in  the  numerical  example  presented 
in  the  Appendix  of  [Blaha,  1989],  illustrating  convergence  properties  of  an 
adjustment  of  a  third-order  polynomial  in  four  variables.  Although  the  standard 
method  converged  slowly  in  one  of  four  analyzed  cases  and  diverged  in  two 
others,  the  geometrical  method  converged  in  two  and  three  iterations, 
respectively.  The  extended  geometrical  method  further  reduced  the  number  of 
iterations  from  three  to  two.  It  is  expected  that  in  most  nonlinear  cases  the 
presence  of  second-order  partial  derivatives  will  translate  into  two  iterations 
in  the  geometrical  method  as  compared  to  several  iterations  needed  by  the 
standard  method. 
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3.  ADJUSTMENT  NITH  NONLINEAR  CONSTRAINTS 


3 . 1  Initial  Relations  in  Matrix  Notatloti 

This  chapter  is  concerned  with  the  introduction  of  nonlinear  constraints 
into  a  nonlinear  least-squares  adjustnent  aodel.  In  standard  adjustment 
notation,  a  set  of  s  constraints  among  u  parameters  is  symbolized  by 

6(X®)  =  0  , 

where  X®  represents  the  set  (column -vector)  of  adjusted  parameters.  The  usual 
approach  consists  in  expanding  these  constraints  in  the  Taylor  series  using  an 
initial  set  of  parametric  values,  symbolized  by  X*^.  It  then  follows  that 

G(X®)  =W  +CX+...=0, 
c 

where 

=  6(X°)  ,  C  -  OG/3X)^  . 

€1  O 

and  where  X=X  -X  is  a  set  of  parametric  corrections.  The  subscript  "o" 

indicates  that  the  matrix  C  is  evaluated  using  the  elements  of  X*^ .  This  matrix 

has  the  dimensions  sxu,  while  the  column  vectors  W  and  X  contain  s  and  u 

c 

elements,  respectively.  In  a  standard  (linearized)  approach,  the  terns 
symbolized  above  by  dots  are  omitted. 

Here,  as  well  as  in  standard  adjustment  theory,  the  constraints  are 
considered  independent,  in  the  sense  that  the  constraint  matrix  C  has  the  full 
row  rank  s.  However,  the  standard  theory  proceeds  in  general  with  By 

contrast,  the  current  nonlinear  development  will  benefit  from  such  initial 
values  X°  for  which  it  holds  true  that 

G(X°)  =  0  .  (8) 

In  this  case,  only  u-s  of  the  u  values  in  can  be  chosen  independently. 
Although  a  set  X°  for  which  G(X°);^0  would  also  be  acceptable,  the  computations 
which  below  will  lead  to  G(x”)=0  would  eventually  have  to  be  performed  as  well, 
and  the  resulting  formulas  would  be  more  cumbersome  and  less  tractable. 

A  simple  iterative  algorithm  leading  to  G(X^)=0  can  be  presented  as 
follows.  First,  the  matrix  C  is  partitioned  into  [C^  C^}.  where  the  submatrices 
and  have  the  dimensions  sx(u-8)  and  sxs,  respectively.  The  submatrix 
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can  be  considered  nonsingular  without  any  loss  of  generality,  since,  if  needed, 

the  parameters  could  always  be  rearranged  beforehand  for  this  condition  to  be 

satisfied.  The  vector  X°  is  similarly  partitioned  into  the  subsets  denoted 

temporarily  and  Xg.  which  contain  u-s  and  s  elements,  respectively.  The 

elements  of  X^  are  chosen  Independently,  and  are  thus  held  fixed  throughout. 

The  remaining  elements,  grouped  in  X^,  are  subject  to  change.  We  symbolize 

their  initial  choice  by  X^  to  accommodate  the  iterative  indices  below.  The 

^  1 

first  iteration  yields  corrections  grouped  in  AX  ,  resulting  in  an  Improved 

1  ^ 
vector  denoted  X^.  After  an  1-th  iteration,  this  vector  becomes 


X 


1 

2 


+  AX 


1 

2  ■ 


If  the  corrections  become  negligible  and  the  iterative  process  is  terminated 
after  n  iterations,  the  values  in  x”  are  adopted  as  the  elements  of  X^-  In 
joining  this  vector  to  the  Independently  chosen  vector  X^ ,  one  obtains  the 
desired  vector  X°. 


3.2  Parawetrlc  Ellalnatlon  Due  to  Constraints 

As  we  have  seen  In  (la,b)  of  Chapter  2,  the  nonlinear  parametric  adjustment 
model  can  be  written  as 

r.a.  r,  a.  .r.a 

X  (u  )  *  X  (u^)  +  A  Au  +  (1/2)  0  *  Au  Au^  +  .  .  .  (9a) 

o  a  (Xp 

where 

C[  Ct  CL 

Au“  =  u"  -  u“  .  (9b) 

with  r  =  l,2 . n.  a  =  l,2 . u;  and  where 

A  =  (3x  /3u  )  n  -  =  (3  X  /3u  3u  )  . 

a  o  <ip  o 

All  lower-case  Greek  letters  vary  in  the  fashion  prescribed  above  for  a.  A 
similar  convention  applies  for  other  kinds  of  indices  as  well  (lower-case  Roman 
letters,  etc.).  The  subscript  ”o”  indicates  the  evaluation  at  the  initial  point 
P  lying  in  the  model  surface. 

The  8  nonlinear  constraints  joined  to  this  model  are  represented  by 

0^(0*)  -  0  .  (10a) 

where  L^^l ,2, . . . ,s.  This  equation  can  be  regarded  as  the  functional  form  of  a 
surface,  generalized  to  higher  dimensions.  In  referring  to  (8)  in  the  preceding 
section,  one  also  has 

G^(u®)  =  0  ,  (10b) 

where  the  values  u",  a=l,2 . u,  are  known.  The  functional  form  (10a) 

o 

restricts  the  final  least-squares  point,  whose  model -surface  coordinates  are  u“, 
to  a  certain  lower-dimensional  surface  embedded  in  the  model  surface.  From 
(10b)  it  follows  that  the  initial  point  P  also  belongs  to  this  lower-dimensional 
surface,  which  will  be  called  "model  subsurface".  This  is  apparent  from  the 
right-hand  sides  of  (10a, b),  which  contain  the  same  sets  of  constants  (zeros). 

The  current  development  is  organized  along  the  following  lines.  First,  the 

coordinate  set  (u*).  a=l,2 . u,  is  partitioned  into  {u^,u^),  A=*l,2 . u-s, 

K°l,2 . s.  This  allows  (10a)  to  be  written  as 


1  2 


0^(0^^,  u*^)  =  0  , 


(11) 


representing  the  functional  forn)  of  tiie  model  subsurface  Subject  to  the 
condition  stated  «;xplici(ly  in  the  sequel,  (11)  makers  it  possible  t()  express  the 
last  s  coordinates  in  terms  of  the  first  u  s  coordinates: 

K  K,  A, 

u  u  (II  )  .  (12) 

Equation  (12)  is  the  Monge  form  of  tlie  model  subsurface  embedded  in  the  model 

surface,  where  u^,  A=1.2 . u-s,  are  the  subsurface  coordinatf:s  (independent 

variables).  The  substitution  of  (12)  into  (11)  yields 

g^(u^)  =  n^(u^,  u'^iu'') )  =  0  .  (13) 

which  is  an  identity  in  the  model  subsurface.  Thus,  further  identities  follow: 

0  ,  3  ^g*' ' =0 .  (14a, b) 

K  A  K 

leading  to  a  relation  for  u  (u  ).  In  this  way,  the  parameters  u  ,  K=1.2 . s, 

will  have  been  effectively  eliminated. 


Expressed  in  the  Taylor  series,  (10a)  reads 


G  ( u  )  5  C  All 
a 


(1'2)  Au“  Au^  ‘ 

ap 


--  0  , 


where  advantage  has  been  taken  of  (lOb),  and  where 


.L 


(3G^'/3u“) 


t-I.  I  ^  2„I. ,  -  o  ^  fi. 

H  „  (3  G  /Oil  3 II  ) 

«P  0 


In  using  the  partition  of  (u  )  and  the  symmetry  of  partial  derivatives  in  the 
lower-case  Greek  Indices,  one  develops  this  equation  as 


A  K,  .  A  ,  K 

ti  ( u  ,  u  )  ?  C .  All  f  C  All 
A  K 


(1.2)  Aii^  Au^ 


o  A  A  K  liL  A  K  ^  M, 

+  2  H.„  All  Au  I  Au  Au  ] 

AK  I\M 


where 


A  ^  A 

Au  -  u 

o 


All 


o 


(15a) 


(15b,cl 


and  where  ii^.  ii*^,  the  model  surface  coordinates  of  P.  are  known  (see  the 
o  o 

relation  10b  and  the  statement  below  it).  Equation  (15a)  corresponds  to  the 
step  represented  by  (11).  In  evoking  (12),  we  next  formulate  the  Taylor  series 
for  II  : 

A  K  A • K  A  A  ,  ,  /o  >  o.K  A  A  .  0 

Au  A  .  Au  ♦  1/2)  0  Au  Au  *  ...  . 

A  A(7 


(16) 


where 


, .  2  K  ^  A  .  Q, 
(o  u  /  dll  dll  ) 


The 


,  ,K  , ^  K , ^  A. 

A '  .  =-  O  u  /  3  U  ) 

A  o 


Q' 


K 


AO  '  'o  '  ■ 

partial  derivatives  are  again  symmetric  in  the  perMneni 


The  substitution  of  (Ifi)  inlo  (Ifia)  in  view  of  the  step 

g'(u‘'')  E  Au'^  •  cj,  A'^  Au*'  -  (1/2)  V.'^  Ai/  ,Mi 

-  (1/2)  Au^  Au^^  -  llj^'  Au^  A'U  Au^' 

A12  AK  w 


J.  ,,K  .A  ..M  .  0 


*  (1/2)  A  '^  All"  All 


0  . 


i nd 1 ces  . 
(13)  yields 

a 


This  identity  is  immediately  confirmed  at  P.  I'pon  differentiating  it  in 
succession  in  accordance  with  the  step  (I4a,b),  and  rearranging  the  free  as  well 
as  the  dummy  indices,  it  follows  that 


2  L,-  A 
3g  /au 


'A 

.1. 


.1.  Q  J,  .  0 

A  ‘  AO  '  ^AO 


,  K  0 
I  A'^  All 


uL  A  V.X 

'"a 


,,l  ,,K  ^  0  ^.M 

»KM  ''  0  A 


(17) 


,2  L  ,  A-  0 
a  g  ,0a  3u 


.  A'  n.X  . 

■  ^ K  ^  AO  ”a0 

«.  H^‘  A  '  ^  A  ' 

KM  0  A 


h''  A'*^  +  A'*^ 

AK  0  OK  A 


(18) 


0 


»(  0 

where  the  dots  represent  terms  containing  Au" ,  Au  Au  .  etc  There  is  no  need  to 
present  partial  derivatives  of  higher  order  than  those  featured  in  (17),  (18). 


The  evaluation  of  (17)  and  (18)  at  P  yields,  respectively. 


where 


MN 


A 


\ 

0 


) 


A ' 


d'*'  pL 

*  L  K 


(19) 

(20) 


The  last  equation  represents  tlie  condition  mentioned  below  (11).  In  matrix 

.M  L 

notation,  this  condition  states  that  the  matrix  (D‘  )  is  the  inverse  of  [C  ]. 

II'  " 

which  In  turn  implies  that  the  matrix  fi;  P 'J  of  dimensions  s^u  roust  have 

A  K 

the  full  row  ranit  s,  and.  therefore,  that  the  constraints  must  be  linearly 


independml.  In  lht>  df'f  i  i-mnl  iv»; ,  nn  evi*i\tiMl  reai'faiif^  j  nji  nf  pat'ameters  will 
pnsura  tliat  the  matrix  I'PHaJar  (i.c..  scjiiarn  and  nonsingular).  This 

subject:  has  already  been  discussed  in  Section  :j.l,  and  has  led  to  (lOb).  Upon 
substituting  (19),  f-'O),  and  liigher  order  partial  derivatives  (not  listed)  into 
(16).  one  obtains  a  relationship  for  u  as  lias  been  indicated  below  (14a. b). 

Next,  An  from  (16)  is  substituted  itito  (9a).  I'pon  ttie  realization  that 

.r  a  .t-.K 

A  Au  A.  Au  +  A,,  Au 
a  A  K 


«  .  „r  ,  A  .  0 


r  .  A  .  K  j- 


K  .  M 


0^-  Au  Au'  n..,  Au‘  Au  ^  2  Au  An  -  Au  Au 


this  substitution  yields 


r  r 
X  .X 

o 


(A^  •  ''^K  ^’a’ 


o  r/’  A  ^  A  A  .  A  A  .  0 

^  V  ^  ()  ••K.M  A  o’ 


where  A''][.  0''^,, . are  kn 

A  AO 
r 

in  (21)  are  interprf‘tpd  as 


.  are  known  from  (19).  (20) 


a:: 


AO  K  AO 


The  symbols  x'  and 


x^(u*)  =  •x'' (  u^  ,  u’^  ( if' )  )  :  x'^(u^) 


r  I',  a,  t  ,  A  K.  A.  -r.  A, 

X  ^  X  ( U  =  X  (  U  ,  U  ( 11  )  )  X  ( u  )  , 

o  o  0  o  o 

p 

where  x  is  known, 
o 

Equation  (21)  is  now  reformulated  to  read 

x''(u^)  x'(u^)  I-  A*  Au^  •  (1/2)  Au'  Au*^  '  .. 

o  A  AO 

where  the  sets  of  impl  i  (M  t  p<j r t  i  a )  ri(* r i  I  i  ves  at  ,  iifi inc'i  y 


(  dx*^,  du^l 


,  3  2  '  r  .  /\  0 

(T  ,x  /dll  .^11  ) 


follow  readily  from  (21): 


a"  A-'^ 
K  ^  A 


(23a) 


"m,  ’  "a„ 


a' 

K  AO 


^  '4  -S 


^KM  ^  A  ''  Q  • 


(23b) 


Whereas  equal irm  (9a)  ifpreseuts  a  nonlinear  model  in  the  pai ameters  u  . 

0  1.2.  .  .n.  e(|U.i'’ion  (22)  represents  .(  nonlinear'  model  in  the  parameters  u'' . 


A  1.2 . II  s  .  In  geomelri  I  <1 1  tt-i-ms,  u  ni-  the  imulel  siirf'iicc  ;  odid  i  nat  ns  nf 

Mu‘  point  depicting  the  unconstrained  least  squares  solution,  and  u  are  the 
model  subsurface  coordinates  of  the  point  depicting  the  constrained  least- 
squares  solution.  The  constiained  solution  can  In-  rallied  out  using  the 
geometrical  algorithms  described  in  Section  2.2  (see  the  geometrical  metliod  or 
the  extended  geometrical  methoti).  Howc-ver.  having  eliminated  the  s  paranu'ters 
u  .  one  now  has  a  smaller  system  to  rt'solve. 


(1  3  Linear  Constraints  as  a  Sper.ial  f.ase 

In  this  section,  we  consider  linear  ronstraints  in  ('onjunct ion  with  the 
parametric  adjustment  model,  whicli  lan  be  citlu'r  linear,  or  nonlinear  as  in 
Chapter  2.  Linear  constraints  (lOa)  would  imply 

,1. 


H 


0 


from  which  it  would  follow  that 
.K 


!!' 

A(l 


0  , 


(24) 


This  outcome  would  lead  to  a  simplification  in  the  formula  (23ij).  where  the 
second  term  in  the  expression  giving  wouid  be  zero, 

Should  the  parametric  model  itself  be  linear,  we  would  further  have 

Q*'-  0 . 

a? 

This,  in  conjunction  with  (24).  would  yield 


"mi  =  " 


(25) 


In  such  a  case,  (22)  would  tun  nnw- 
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Fqual ions  (26a. b)  represent  a  linear  parametrii  adjustment  model,  where  the 

original  design  tensor  a'  is  replaced  l>y  and  the  original  set  An”  is 
A  ' 

replaced  by  Au  .  As  in  the  preceding  sect i mi.  this  system  is  smaller  than  the 
origitiiil  one.  diiP  to  the  elimination  of  the  parameteis  n 


The  above  eliminatioti  eaii  be  eoiil  irmed  in  tfie  standard  adjustment  notation 


as  follows.  In  a  linear  mo<lel ,  the  exact  observation  equations  read 

V  -  A  X  ‘  r,  .  (27) 

where  I,-L^  and  X^x’-x”.  as  defined  earlier.  The  design  matrix  A  is  now 
constant  regardless  of  the  set  X^'  The  linear  constraints  are  expressed  by 

G(X^)  =-  p  ‘  r  X'”'  -  0  .  (28a) 

where  p  is  a  known  constant  set  <jf  s  elements  and.  in  analogy  to  A,  the  matrix  C 
is  constant  Iti  agreement  with  (8).  we  us<*  tht,‘  ijiitial  values  x'^  siu-h  that 

G(x'’)  -  p  -  r  x“  0  .  (2Hb) 

and  state  ttiat  only  a  us  subset  of  x”  can  be  chosen  independently  in 
partitioning  r  as  in  Section  3.1.  i.e  ,  r:=(f'^  partitioning  X°  similarly 

as  X^5[x”^  where  X*^  is  the  chos<>ti  subset,  from  (28b)  we  have 

x"  ^ 


Due  to  (28a.b),  we  can  write 
r  X  --  0  . 


(30) 


In  partitioning  X  fx|^  X^)^  from  (30)  we  deduce,  similar  to  (29): 

x,  =  r,'  c,  X, 


(31  ) 


If  A  in  (27)  is  partitione<l  in  accordance  with  X  as  [A^  ^2^’  follows  that 


v. 


V  -  A  X  ‘  A  X  ‘  L 
1  1  2  2 


However,  the  substitution  of  (31)  into  this  relation  yields 


f.  . 


(32a) 


where 


1 


N  ^2  ^^2  S 


(32b) 


The  matrices  and  have  already  been  used  in  forming  in  (29).  We 

observe  that  the  system  (3,2a,b)  corresponds  to  (2Gb)  together  with  (19). 


4.  CONCLUSION 


Many  geodetic  problems  aie  either  presented  in  the  form  of  a  parametric 
model,  or  can  acquire  this  form  upon  a  simple  linear  transformation  of 
observables.  A  good  example  of  this  transi orraat ion  is  offered  by  the  SST 
adjustment  model  The  variance  covariance  matrix  of  the  original  observations 
in  such  cases  must  be  transform(?d  accordingly,  if  the  rigor  of  the  adjustment  is 
not  to  be  compromised,  The  most  widely  a< i epicd  method  of  afijustment,  the 
least -squares  method,  is  used  in  piactice  for  motiels  that  are  linear,  or  have 
been  linearized.  By  contrast,  the  least  -  squares  approach  presented  herein 
foetuses  on  the  nonlinear  parametric  adjustment  model  which  may,  furthermore, 
contain  a  set  of  nonlinear  constraints  among  the  parameters. 

The  resolution  of  the  nonlinear  parametric  adjustment  model  without 
constraints  is  addressed  through  an  isomorphic  geometrical  setup  with  tensor 
structure  and  notation,  represented  by  a  u  dimensional  model  surface  embedded  in 
a  flat  n-dimensional  observational  space.  The  n  observations  correspond  to  the 
observational -space  coordinates  of  the  point  Q,  the  u  initial  parameters 
correspond  to  the  model-surface  coordinates  of  the  initial  point  P,  and  the  u 
adjusted  parameters  correspond  to  the  model  surface  coordinates  of  the  least - 
squares  point  P.  The  least  squares  criterion  results  in  a  minimum  distance 
property  implying  that  the  vector  PQ  must  be  oithogonal  to  the  model  surface. 

The  geometrical  setup  leads  to  the  solution  of  modified  normal  equations, 
characterized  by  a  positive  definite  matrix  The  latter  cont.ains  second  order 
and.  optionally,  third-order  partial  derivatives  of  the  observables  with  respect 
to  the  parameters.  This  approach  significantly  shortens  the  convergence  process 
as  compared  to  the  stajidard  (linearized)  method. 

The  nonlinear  paranieti ic  adjustment  model  with  nonlinear  constraints  is 
also  resolved  through  geometrical  analogies  In  this  situation,  a  point 
representing  the  least  squares  solution  is  restiict<;d  to  lie  in  the  model 
subsurface,  i.e.,  a  surface  of  smaller  dimeusions  than  the  model  surface  in 
which  it  is  embedded.  The  geometrical  approach  leads  to  the  replacement  of  the 
model  surface  by  the  model  subsurface,  and  to  the  treatment  of  the  observational 
point  Q  with  respect  to  the  new  stuface  in  tin-  manner  that  resulted  in  the 
( unrest  r  i  cted )  least  s(piai-es  point  P  Ai  cord  ingl  y ,  the  constrained  least 
squares  point  is  the  result  of  an  orthogonal  prejeeijou  of  g  onto  the  model 


subsurfacR.  In  the  adjusVnnMit  i iiology .  this  approach  eliminat  s  s  of  the 

original  u  parameters,  wliere  s  is  alsd  the  numbei  of  constraints.  The  temaining 
parameters  are  resolveil  by  the  method  of  the  nonlinear  parametric  least  squart-s 
adjustmc'nt  without  (constraints,  wiier*'  all  the  arrays  must  be  proptrrJy  m  '  lificd. 

A  special  (class  of  constCftints  is  represented  by  inequality  constraints, 
which  art!  treated  here  in  a  lintcai  form,  and  arc  cousidered  in  conjunction  with 
a  linear  parametric  adJrrslmiMit  model  The  isomorjili  i  c  (p'ornetrical  setup  is  now 
partly  simplifieil,  in  the  stcnse  that  geiier.il  surfarc-.  .rrc  replaced  by 
hyperplanes.  The  topi(  of  rneipialirv  constr.iints  diffei's  from  its  e(|uali(y 
counttrrpart  in  that  only  some  constraints  (terivied  binding)  arte  i-etained  and 
subsequently  enfoin  ('d  as  equ.iltty  constraints,  wtreicas  the  remaining  constraints 
(called  nonh  i  nd  ing )  arte  ignored.  The  most  difficult  (ju(‘st  i  on .  then,  is  to 
determine  wltich  of  the  constraints  siioiiid  be  retained  as  binding.  In  the 
geometrical  context,  this  pidhlcm  is  addressed  by  orthogonally  projecting  the 
point  reprrrsent  i  ng  the  urrrest  r  i  I't  ed  least  srfuares  soitttion  from  an  origirral 
model  hyperplane  onto  appropriate  hyperplancs  of  lower'  dimertsions  Since 
orthogon.il  projections  resirlt  in  the  shortest  possible  distance  between  the 
unrestricted  least  .sqttares  poirrt  and  the  final  corrst  ra  ined  point,  the  solutiorr 
belongs  to  the  least  sqrrares  category 

The  geomeli  i  c-r  1  approach  to  i  rterpi.r  1  i  ty  cotts  t  ra  i  ii  t  s  is  (compared  with  the 
standard  resolution  of  the  same  least  s(iuares  problem  via  quadratic  programming. 
A  nrrmerical  e.x.tmple  with  four  p.rtameteis  is  solved  by  both  methods,  leading  to 
identical  results.  The  main  diffecenre  betwt'en  the  two  methods  is  conceptual. 
The  former  attributes  a  clearsnt  geortic  t  r  t  ca  I  meaning  to  every  ad.iirstment 
quantify,  and  re.o  hes  the  rorrst  r.i  itted  le.tst  squares  solirtion  in  accordatt(ce  with 
geometrical  principles  The  lalft-r  is  based  on  algebraic  principles,  the  key 
element  of  which  is  the  Gauss  .lord.iti  pivoting,  (hi  f  h('  oper.afional  level,  the 
numerical  systems  treated  by  the  geometrical  algorithm  become  progressively 
smaller  after  each  orthogonal  projection  along  a  given  patli.  By  contrast,  the 
size  of  such  systems  treated  by  the  quadratic  programming  algorithm  rem.iins 
constant.  Another  advantage  of  the  geometriial  ajrproach  consists  in  an  early 
detection  of  a  path  leading  to  a  non  sriuares  solution. 


APPENDIX  A 


RIGOROUS  ADJUSTMENT 

OF  SATELLITE-TO-SATBLLITE  TRACKING  DATA 

In  this  appendix,  we  describe  the  adjustaent  aodel  of  the  satelllte-to- 
satelllte  tracking  (SST),  and.  subsequently,  a  nonlinear  least-squares  Method 
consistent  with  such  a  model.  The  observations  in  the  SST.  whether  in  the  high- 
low  or  the  low-low  configurations,  are  the  relative  velocities  between  two 
satellites.  These  range-rate  data,  the  error  characteristics  of  which  are 
assumed  to  be  known,  serve  in  the  deteraination  of  the  detailed  gravity  field  of 
the  earth. 


SST  Adlustaent  Model 

In  terms  of  the  adjustment  model,  the  observables  are  represented  by 
Intersatellite  range  rates  and  the  unknown  parameters  are  represented  by 
selected  gravity  field  parameters  and  other  desired  quantities.  In  a  standard 
procedure,  applied  to  a  variety  of  problems  in  physical  sciences,  the  adjustment 
model  Is  linearized,  whereupon  range-rate  observations  give  rise  to  observation 
equations.  This  system  can  then  be  resolved  by  the  parametric  method  (also 
called  the  observation  equation  method)  of  the  least-squares  theory.  It  should 
be  pointed  out,  however,  that  the  intersalel lite  range  rate  does  not  provide  a 
direct  measure  of  the  potential  at  the  satf'lllte  positions.  A  more  direct 
relationship  would  be  obtained  if  the  observables  were  Intersatell Ite  velocity 
rates,  i.e.,  relative  accelerations  between  two  satellites. 

This  point  is  illustrated  in  inertial  Cartesian  coordinates  with  familiar 
vector  notations  (here  the  vectors  are  underlined).  With  1=1,2,  is  the 
position  vector  of  the  satellite  1,  is  its  velocity  vector.  R  is  the 
magnitude  of  the  relative-position  vector  between  the  two  satellites,  e  is 
their  unit  relative-position  vector,  and  R  is  their  range  rate.  The  latter  is 
the  projection  of  the  relative  velocity  vector  onto  the  relative  position 
vector : 


-  Xj )*e  . 

e  -  (X^ 

-  Xj)/R  . 

R  =  IX^ 
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The  quantity  central  to  the  deterninatlon  of  the  gravity  field  is  the  time 
derivative  of  R,  obtained  by  numerical  differentiation.  The  new  quantity  R, 
i.e..  the  relative  acceleration  between  the  two  satellites,  can  be  directly 
related  to  the  potential  at  satellite  positions.  The  acceleration  vector  is 
numerically  equal  to  the  gradient  of  the  gravitational  potential.  We  thus  have 

R  2  dR/dt  =  a  V  b  , 

where 

a  -  -  ^j)-e  . 

b  =  [(Xg  '  -  Xj)  -  R^l/R  . 

Here  denotes  the  potential  at  the  satellite  position  1  and  is  the 
gradient  of  .  The  main  contribution  of  the  relative  acceleration  between  the 
two  satellites  is  contained  in  the  term  a,  which  is  the  projection  of  the 
relative-acceleration  vector  onto  the  relative-position  vector. 

Although  the  relative  acceleration  R  is  highly  suitable  for  modeling  the 
gravitational  potential  as  evidenced  by  the  term  a  above,  it  is  not  an  observed 
but  a  derived  quantity.  It  Is  usually  obtained  by  a  curve-fitting  procedure, 
i.e.,  by  the  filtering  and  the  smoothing  of  the  original  observations  R, 
followed  by  a  numerical  differentiation  with  respect  to  time.  There  exist 
different  methods,  such  as  the  spline  method,  for  converting  the  original 
relative-velocity  observations  R  into  the  relative-acceleration  "observations" 
R.  Various  aspects  of  such  methods  are  described,  for  example,  in  [Rummel  et 
al . ,  1976]  and  in  [Hajela,  1977].  During  this  mathematical  treatment,  the 
original  data  are  modified  (by  curve  fitting),  and  the  modified  data  are  then 
transformed  Into  quantities  of  a  different  kind  (by  numerical  differentiation). 
The  resulting  data  are  not  unique,  i.e.,  a  different  data  set  is  obtained  with 
each  different  method  used.  As  a  consequence  of  the  accumulated  modifications 
of  the  original  data,  the  original  error  characteristics  are  lost.  The 
variance-covariance  matrix  of  the  new  data  set  entering  the  least-squares 
adjustment  must  be  then  supplied  in  some  approximate  fashion. 

A  more  rigorous  approach  would  be  to  utilize  the  original  observations 
together  with  their  variance-covariance  matrix.  This  can  be  done  by 
differencing  the  range-rate  observations  at  some  suitable  time  intervals.  In 
particular,  one  has 


21 


[R(t*^At)  -  R(t)]/At  =  R(t)  +  O  , 
o  -  (l/2)(dR/dt) At  V  . .  .  . 


In  considering  the  previous  outcone,  the  SST  adjustment  model  becomes 

[R(t  +  At)  -  R(t)]/At  =  (^2  ■  +  o  +  b  .  (1) 

This  model  relates  a  combination  of  observations  (two  per  equation)  to  the 
parameters  of  the  gravity  field  and  to  the  small  corrections  o  and  b. 


Transformation  of  the  SST  Adjustment  Model 

An  adjustment  model  represented  by  a  system  of  equations  encompassing  both 
the  observables  and  the  parameters  is  described  in  general  by 

f(X®.L®)  =  0  ,  (2) 

d  O  d  b 

where,  according  to  the  standard  adjustment  notation,  X  =X  ♦X  and  L  -L  >V.  The 

d  o 

sets  (column-vectors)  X  ,  X  ,  and  X  contain  the  values  of  adjusted  parameters, 
initial  parameters,  and  parametric  corrections,  respectively;  and  the  sets  L  , 
L^,  and  V  represent  adjusted  observations,  actual  observations,  and  residuals, 
respectively.  It  is  an  ongoing  practice  that  a  model  such  as  (2)  is  linearized 
upon  neglecting  higher-order  terms  in  its  Taylor  series  expansion.  This  gives 
rise  to  the  matrix  equation 

AX+BV+W=0,  (3) 

where  A  =  3f/0X  and  B  =  0f/aL,  both  evaluated  with  X°  and  L^,  and  where 
W=f(X°,L^).  The  matrices  A  and  B  are  assumed  to  have  the  full  column  rank  and 
the  full  row  rank,  respectively.  Equation  (3)  characterizes  the  standard  setup 
of  the  general  adjustment  method,  which  is  resolved  in  accordance  with  the 
least-squares  principle. 

The  neglect  of  higher-order  terms  In  the  Taylor  series  expansion  of  a 
nonlinear  model  represents  the  greatest  simplification  and,  at  the  same  tine, 
the  greatest  shortcoming  of  the  standard  adjustment  theory.  The  price  to  pay 
for  such  a  simplification  Is  represented  either  by  non-rigorous  results  if  the 
solution  is  not  iterated,  or  by  the  necessity  to  iterate  the  least-squares 
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algorithm.  Depending  on  how  severe  is  the  model's  nonlinearity,  the  Iterative 
process  may  be  slow  to  converge.  Thus,  in  theory,  equation  (3)  should  read 

A,XfBV  +  W+  ...=0,  (4) 

where  the  dots  represent  the  contribution  of  nonlinear  terms  due,  in  general,  to 
the  model's  nonlinearity  in  both  the  observables  and  the  parameters. 

If,  however,  the  observables  (but  not  the  parameters)  are  combined  in  a 
linear  fashion,  (2)  can  be  written  in  a  different  functional  form,  namely 

-  B  =  F(X''')  .  (5) 

“  3 

where  the  symbol  L  represents  "transformed  adjusted  observations".  The  new 

3  b 

matrix  B  is  again  assumed  to  have  the  full  row  rank.  Since  L  =L  +V  and 
L^=L*^^V,  where  symbolizes  "transformed  observations"  and  V  symbolizes 

"transformed  residuals",  one  has 

=  B  ,  V  =  B  V  .  (6a, b) 

If  the  model  (5)  is  now  linearized  in  the  parameters,  it  results  in  the 
following  system  of  observation  equations: 

V  =  A  X  L  ,  (7) 

where 

A  =  3F/9X  =  design  matrix  (evaluated  with  the  initial  set  X°), 

X  =  X®  -  X°  =  vector  of  parametric  corrections,  and 
L  =  F(X°)  -  i,^  =  F(X°)  -  B  =  vector  of  constant  terms. 

The  variance-covariance  matrix  £  of  transformed  observations  is 
formed  rigorously  from  (6a)  as 

£  =  B  £  ,  (8) 

where  £  is  the  variance-covariance  matrix  of  the  actual  observations  L^.  The 
adjustment  of  the  linearized  parametric  model  (7)  would  proceed  with  the  proper 
variance-covariance  matrix  from  (8),  whose  inverse  would  be  the  weight  matrix  of 
transformed  observations  and  would  be  used  in  the  formation  of  normal  equations. 
The  shortcoming  associated  with  the  linearization  of  a  nonlinear  adjustment 
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model  has  been  discussed  at  the  outset  of  this  section.  Thus,  In  theory, 
equation  (7)  should  read 

V=AX+...+L.  (9) 

where  the  dots  represent  the  contribution  of  nonlinear  terms  due  to  the  model's 
nonlinearity  in  the  parameters. 


The  SST  adjustment  model  presented  in  (1)  has  the  form  equivalent  to  (9). 
This  stems  from  the  fact  that  the  observations  represented  by  the  left-hand  side 
of  (1)  conform  to  the  linear  pattern  (6a).  In  particular,  the  left-hand  side  of 
an  i-th  equation  is  formed  as  ( 1 /At  )x  (observat  ion  1*^1)  ( 1/At)x  (observation  i). 

With  a  constant  At,  the  matrix  B  would  have  the  form 


B 


(1/At) 


-1  1 

0  1 


0  0 
1  0 


With  the  aid  of  this  matrix  and  of  the  rigorous  variance  covariance  matrix  I  of 
range -rate  measurements,  one  forms  the  rigorous  variance-covariance  matrix  £ 
as  indicated  in  (8),  which  is  then  to  be  used  in  the  parametric  least -squares 
adjustment  of  the  SST  model. 

The  SST  model  is  in  general  nonlinear  in  the  parameters,  with  the  degree  of 
nonlinearity  depending  on  the  type  of  parameters  expressing  the  desired 
components  of  the  earth's  gravity  field  and  other  phenomena.  This,  together 
with  the  above  outcome,  has  compelled  the  analysis  to  shift  focus  from  the 
general  adjustment  model  to  the  nonlinear  parametric  adjustment  model.  With  the 
provision  of  using  the  proper  variance-covariance  matrix  from  (8),  the  overbars 
are  dropped  and  the  parametric  model  is  written  as 

L®  =  F(X^)  =  F(X°)  ^  A  X  ♦  . . .  (10) 

This  model  is  alternately  presented  in  the  form  (9),  which  now  reads 

V  =  A  X  .  (11 ) 

where 

V  =  I,‘*  -  =  vector  of  residuals, 

A  =  3F/0X  =  design  matrix  (evaluated  with  the  initial  set  ) , 
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X  =  ^  vector  of  parametric  corrections,  and 

L  =  F{X°)  -  ="  vector  of  constant  terms. 

The  variance-covariance  and  the  weight  matrices  of  observations  are  denoted  Z 

*  1  £1 
and  P,  respectively,  where  P=Z  There  exist  an  infinite  number  of  sets  L 

consistent  with  the  model  (10).  Of  these,  the  least-squares  principle  selects 

the  one  fulfilling  V^PV^^minlmum. 


->5 


APPENDIX  B 


LINEAR  PARAMETRIC  ADJUSTMENT  MODEL 
MITH  LINEAR  INEQUALITY  CONSTRAINTS 

Introduction  to  Adjustment  with  Ineqiutiity  Constraints 

When  used  in  conjunction  with  a  parametric  adjustment  model,  inequality 
constraints  limit  the  domain  of  selected  parameters  or  functions  thereof. 
Qualitatively,  this  property  reminds  one  of  the  familiar  equality  constraints, 
where  the  individual  signs  -  would  now  be  replaced  by  the  signs  )  or  < . 

However,  a  set  of  inequality  constraints  is  in  general  less  restrictive  than  its 
equality  counterpart  because  only  a  subset  of  the  former  changes  into  a  subset 
of  the  latter  in  the  course  of  adjustment  computations.  The  constraints  having 
this  property  (i.e..  whose  signs  >  or  <  have  changed  into  the  sign  =)  are  called 
binding:  their  number  may  range  from  zero  to  the  total  number  of  the  original 
constraints.  The  remaining  constraints  are  called  nonbinding,  and  they  have  no 
effect  on  the  adjustment. 

The  basic,  and  perhaps  the  most  common,  Itind  of  inequality  constraints 
consists  of  upper  or  lower  bounds  Imposed  on  linear  combinations  of  parameters, 
including  the  class  of  upper  or  lower  bounds  imposed  on  selected  parameters 
themselves.  An  Imposition  of  known  lower  bounds  on  selected  parameters,  for 
example,  can  be  an  Important  asset  if  this  is  dictated  by  physical  or 
mathematical  reality.  One  cannot  Include  such  vague  information  into  the 
adjustment  process  by  weighting  the  parameters  in  question,  since  this  would 
necessitate  their  a  priori  estimates  as  well  as  the  variance-covariance  matrix 
associated  with  such  estimates.  On  the  other  hand,  discarding  this  or  similar 
information  might  lead  to  results  Inconsistent  with  tiie  reality.  One  can  then 
decide  to  either  forego  the  inclusion  of  inequality  constraints  in  the  hope  that 
such  an  inconsistency  will  be  small  or  nonexistent  due  to  appropriate  measuring 
and  modeling  techniques,  or  to  incorporate  these  constraints  into  a  rigorous 
least -squares  adjustment  and  thereby  transform  the  hope  Into  certainty. 

Adjustment  using  bounds  on  specific  parameters  or  their  combinations  Is 
Increasingly  finding  Its  niche  in  geodesy,  photogrammetry .  oceanography,  and  In 
many  other  sciences.  For  example,  one  of  the  most  recent  photogrammetrlc 


applications  is  concerned  with  digital  object  reconstruction.  It  has  been  noted 
that  the  standard  (lineai)  least  squares  method  is  often  insufficient  because  of 
its  smoothing  effect.  A  suitable  alternative  has  been  considered  in  terms  of 
"Chebyshev  formulation",  which  represents  a  generalization  of  this  method.  In 
recent  reports  and  papers,  such  as  [Frltsch,  1987],  this  formulation  is  achieved 
by  using  the  concept  of  least-squares  adjustment  with  inequality  constraints, 
where  both  the  parametric  model  and  the  constraints  are  linear.  The  purpose  of 
such  a  development  is  to  reduce  the  solution  space  to  a  band,  which,  in  turn, 
allows  one  to  minimize  the  maximum  error,  l.e.,  to  accommodate  the  worst  case. 
According  to  Fritsch  [1987],  the  inequality-constrained  least-squares  adjustment 
was  introduced  into  geodesy  by  B.  Schaffrin  in  1981. 

As  a  plausible  oceanographic  illustration  with  bounds  imposed  on  the 
parameters,  we  mention  the  action  spectral  density  of  fluctuations  of  sea 
surface  elevation,  which  is  a  non-negative  quantity  over  all  spectral  bands. 

The  sea  surface  fluctuations  have  been  studied  by  Snyder  [1988],  who  utilizes 
the  linear  relation  A=ZA,G..  where  A  is  the  action  spectral  density,  A,  are 
parameters  (to  be  determined  from  observations),  are  basis  functions,  and 

where  the  summation  extends  over  the  spectral  bands  considered,  i=l,2 .  In 

a  convenient  approach,  the  spectra]  representation  is  adopted  as  piecewise 
continuous,  such  that  Inside  the  i  th  spectral  band  and  G^=0  elsewhere. 

The  key  consideration  pertinent  to  our  discussion  resides  in  the  fact  that 
should  A>0  hold  true  everywhere,  all  of  the  parameters  A^  are  required  to  be 

non-negative.  Thus,  the  condition  A.>0,  i=l,2 .  represents  a  basic  case  of 

linear  inequality  constraints. 

Motivated  by  the  above  considerations,  the  present  study  has  focused  on 
linear  inequality  constraints.  Due  to  complexities  associated  with  nonlinear 
least -squares  adjustment,  certain  aspects  of  which  are  treated  in  [Blaha,  1989], 
the  analysis  has  been  further  restricted  to  applications  involving  linear 
parametric  adjustment.  Upon  taking  advantage  of  a  geometrical  setup  reflecting 
the  situation  where  a  linear  model  is  to  be  resolved  in  conjunction  with  linear 
inequality  constraints,  a  new.  yet  relatively  simple  algorithm  has  been  derived 
producing  a  unique  least  squares  estimate.  A  similar  isomorphic  geometrical 
setup  can  undoubtedly  serve  in  the  future  in  resolving  also  nonlinear  adjustment 
models  in  conjunction  with  linear  and  even  nonlinear  Inequality  constraints. 


Matrix  ForiBulation  of  Linear  InetiuHl  ity  Constraints 


In  this  section,  we  present  soiie  of  the  outcome  of  the  study  which  will  be 
described  in  the  AFGL  Scientific  Report  No  2,  Linear  Paiametric  Adjustment 
Model  with  Linear  Equality  and  Inequal  it  y  constraints.  Although  the  algorithm 
for  resolving  the  linear  adjustment  model  with  linear  inequality  constraints  has 
b(,'en  derived  using  geometry  with  tensor  structure  and  notation,  here  the  results 
are  presented  in  the  standard  matrix  notation  The  inequality  constraints  have 
the  form  of  lower  bounds,  since  upper-bound  constraints  can  be  transformed  into 
lower-bound  constraints  upon  multiplying  the  pertinent  inequality  by  -1. 

The  linear  parametric  adjustment  model,  before  the  introduction  of  any 
constraints,  reads 

_a  ^b  ,, 

L  =  L  +  V=  AXfL  , 

where  the  column-vector  L'*  contains  n  adjusted  observations.  contains  n 
actual  observations,  V  contains  n  residuals,  L*’  contains  n  constant  values,  and 
the  column  vector  X  contains  u  parameters,  u<n  The  symbol  A  denotes  the  design 
matrix  of  dimensions  nxu,  assumed  to  have  the  full  column  rank  u.  The  above 
relation  is  often  written  in  the  form  of  observation  equations 

V  =  A  X  +  I,  ,  (1) 

where  L=b°-L^.  When  subjected  to  the  least  squares  criterion 

P  V  =  minimum  ,  (2) 


where  P  is  the  weight  matrix  of  observations  adopted  as  the  inverse  of  E,  the 
variance-covariance  matrix  of  observations.  (1)  yields  the  normal  equations 


N  X  --  U 


(3a) 


with 


N  -  A^  P  A  ,  (I  -  P  1,  ;  (3b.c) 

here  N  is  a  positive  definite  matrix  of  dimensions  uxu  and  U  is  a  column -vector 
of  u  elements. 


The  scope  of  a  study  conct^rned  with  linear  inequality  constraints  can  be 
narrowed  down  In  two  ways.  First,  the  constraints  can  be  treated  in  the  form 

(  X  3  0  .  (4) 


O  C; 


where  C  is  a  matrix  of  dimensions  s^u,  with  s<u,  assumed  to  liave  the  full  row 
rank  s.  Since,  in  general,  X*^,  where  X^  symbolizes  the  adjusted 

parameters  and  symbolizes  the  initial  values  of  parameters,  the  cases  such  as 
CX>c,  where  e  is  a  constant  vector  of  s  elements,  can  be  transformed  into  (4) 
upon  properly  modifying  the  values  in  X^*. 

The  second  simplification  can  be  achieved  through  a  unique  linear 
parametric  transformation  carrying  tier  vector  X  into  a  vector  Y,  likewise 
composed  of  u  elements.  In  particular,  upon  partitioning  X  and  Y  into  u-s  and  s 
elements,  and  attributing  the  .sut)sets  a  prime  and  a  double  prime,  respectively, 
we  have 


Y'  -  X'  , 

Y”  =  C  X'  f  C"  X"  ^  0  , 

where  C  has  similarly  been  partitioned  into  the  submatrices  C  and  C" .  The 
second  equation  above  is  equivalent  to  (4).  From  these  relations  one  can 
express  X'  and  X"  in  terms  of  Y'  and  Y",  and  use  the  result  in  (1).  yielding 


V  =  A’  Y  ^  I, 


where  A'  follows  from  A  and  from  the?  transformation  coefficients.  We  note  that 
C"  has  been  assumed  to  be  a  regular  matrix  (i.e..  square  and  nonsingular).  Due 
to  the  full  row  rank  of  C,  this  assumption  is  justified,  either  initially  or 
upon  renumbering  the  parameters.  After  the  vector  Y  and  its  variance-covariance 
matrix  have  been  determined,  the  parametric  transformation  yields  the 
original  vector  X,  and  the  law  of  variance-covariance  propagation  yields  . 

In  order  to  reduce  the  number  of  symbols,  we  change  the  notation  from  A'  to 
A,  and  from  Y  to  X.  The  last  two  ecpiations  are  then  transcribed  as 


V  -  A  X  h  .  {5a) 

X"  ^  0  ,  (5b) 


where  X"  is  a  subset  of  X.  If  X"  comprises  the  fuli  set  X,  tlie  adjustment  model 
with  constraints  becomes 


V  -  A  X  ‘  I,  ,  (6a) 

X  ^  n  .  (6b) 


This  Is.  in  fact,  the  basic  pioblem  discussed  in  the  previous  section,  where  all 
the  parameters  were  required  to  be  non-negative.  However,  this  case  is  quite 
suitable  for  a  general  treatment  of  linear  inequality  constraints,  since 
disregarding  the  pertinent  n-s  constraints  in  (6b)  leads  to  the  system  (5a, b), 
which,  in  turn,  is  equivalent  to  the  system  (1),  (4). 


Outline  of  the  Geometrical  Approach 

The  derivations  summarized  in  this  section  will  appear  in  the  above-cited 
report.  In  the  geometrical  context,  the  n  observations  grouped  in  correspond 
to  the  coordinates  of  the  "observational  point"  Q  lying  in  an  n  dimensional  flat 
"observational  space".  Since  the  metric  tensors  of  all  the  manifolds  considered 
here  are  constant  (due  to  the  linear  setup),  we  can  interpret  coordinates  of  a 
point  as  contravariant  components  of  its  position  vector,  and  coordinate  lines 
as  oblique  Cartesian  axes  with  constant  individual  scales  The  matrices  Z  and  P 
(constant)  correspond  respectively  to  the  associated  metric  tensor  and  to  the 
metric  tensor  of  the  observational  space.  The  u  initial  parametric  values  in  X° 
correspond  to  the  coordinates  of  the  "initial  point"  P^.  The  latter  lies  in  the 
known  u-dimensional  flat  model  surface",  also  called  "model  hyperplane",  which 
is  embedded  in  the  observational  space.  In  an  unrestricted  least-squares  (L.S) 
adjustment,  the  adjusted  parameters  X**  correspond  to  the  point  denoted  P,  which 
likewise  lies  in  the  model  surface.  The  adjustment  notation  X  then  designates 
the  contravariant  components  of  the  vector  P^P-  Hue  to  constant  metric  tensors 
of  the  observational  space  and  of  the  model  surface,  vectors  in  these  manifolds 
can  be  freely  parallel-transported  to  any  location.  This  allows  us  to  identify 
P^.  throughout  the  analysis,  with  the  coordinate  origin  in  the  model  surface. 

The  unrestricted  L.S.  solution  represented  by  the  point  P  is  obtained  by 
projecting  the  observational  point  Q  onto  the  known  u  dimensional  model 
hyperplane.  Here  the  term  "projection"  will  always  be  synonymous  with 
"orthogonal  projection".  The  L.S.  solution  subject  to  inequality  constraints 
would  be  obtained  by  projecting  Q  onto  another  surface,  as  yet  unknown,  and 
generating  the  point  denoted  P.  Clearly,  in  the  absence  of  constraints,  or  in 
the  presence  of  only  nonbinding  constraints,  P  and  P  would  coincide.  But 
since,  in  general,  the  constraints  (5b)  or  (6b)  limit  the  solution  point  to  a 
region  of  the  ii  dimensional  model  hyperplane,  the  point  P  must  be  transferred  in 
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so*e  way  Into  this  "admissible"  region.  At  the  same  time,  iji  order  to  fulfill 
the  [,.S.  criterion,  the  new  point  must  coincide  with  P. 

The  above  discussion  indicates  that  the  point  P,  itself  the  result  of  a 
projection,  must  be  further  projected  onto  an  "envelope"  delimiting  the 
admissible  region.  The  transfer  of  P  inside  this  region,  or  onto  its  envelope 
other  than  by  a  projection,  would  be  inconsistent  with  the  L.S.  principle, 
which,  in  the  geometrical  conte,xt,  translates  into  the  shortest  distance 
principle.  The  form  of  (v5b)  or  (fib)  allows  us  to  consider  the  "sides"  of  the 
above  envelope  as  portions  of  lower  dimensional  hyperplanes  embedded  In  the 
model  hyperplane  and  spanned  by  combinations  of  coordinate  axes.  We  will  solve 
the  Inequality-constrained  L.S.  adjustment  by  projecting  the  unrestricted  L.S. 
point  P  from  the  known  u  dimensional  model  hyperplane  onto  a  lower-dimensional 
"embedded  hyperplane",  thereby  generi>ting  the  point  P,  in  such  a  way  that 

a)  the  latter  hyperplane  has  the  highest  dimensionality  possible: 

b)  the  linf;  connecting  P  and  P  does  not  pass  through  any  part  of  the 
admissible  region;  and 

c)  the  point  P  is  consistent  with  the  inequality  constraints. 

The  thus  generated  constrained  L.S.  point  P  is  identical  to  the  point  which 
would  be  obtained  via  an  orthogonal  projection  of  the  observational  point  Q 
directly  onto  the  (unknown)  lower  dimensional  hyperplane. 

The  determination  of  the  unrestricted  L.S.  point  P  corresponds,  in  the 
adjustment  context,  to  the  solution  of  the  normal  equations  (3a-c).  If  the 
parameters  sliow  conflict  with  the  constraints  represented  by  (5b)  or  (fib),  P 
must  be  projected  as  discussed  above.  The  geometrical  algorithm  developed  for 
this  purpose  proceeds  in  accordani;e  with  the  strategy  summarized  as  follows: 

1)  Only  the  conflicting  (negative)  elements  of  X  will  induce  projections. 
(Since  It  can  be  shown  that  any  permutation  of  a  given  sequence  of  projections 
yields  the  same  point,  a  sequence  where  only  negative  components  have  induced 
projections  can  be  imagined  In  a  different  permutation,  which  now  may  include 
projections  corresponding  to  positive  components.  This  specific  sequence  would 
also  yield  the  above  point,  but  sucli  arrangements  are  strictly  avoided.) 

2)  The  projections  are  "closely  nested",  in  th*-  sense  that  a  point  lying  in  a 
given  hyperplane  can  only  be  projected  onto  a  liypetplane  whose  dimensionality  i 


lower  by  one.  This  process  is  repeated  until  the  constrained  L.s.  point  P  Is 
reached,  or  until  a  given  sequence  is  rejected. 

The  point  P  can  be  reached  along  different  routes,  or  "branches",  formed 
by  allowable  permutations  in  a  successful  sequence  of  projections.  Similarly, 
one  or  more  non  b.S.  points  (to  be  rejected)  may  be  reached  by  many  branches. 

It  is  thus  important  to  avoid  branches  that  are  essentially  repetitious.  To 
this  end,  the  geometrical  algorithm  will  keep  track,  at  each  level,  of  the 
permutations  that  have  already  taken  place.  The  level  m  is  defined  as  the  stage 
in  the  algorithm  where  the  number  of  consecutive  projections  has  reached  m,  and 
the  dimensionality  of  the  latest  hyperplane  is  urn.  If,  for  example,  a  given 
branch  includes  the  projections  symbolized  by  the  letters  a,  b,  and  c,  another 
branch  containing  the  projections  c,  a,  b  will  be  discarded  at  the  third  level. 
As  anothei  example,  if  the  projection  sequence  g.  h  is  rejected,  the  sequence  h. 
p.  q.  g  will  be  automatically  rejected  since  it  contains  the  forbidden 
combination  g  and  h. 

One  of  the  crucial  elements  in  the  geometrical  algorithm  is  the  early 
detection  of  branches  that  should  be  rejected.  The  rejection  criterion  will  be 
provided  by  a  "guide  vector"  of  u  elements,  which  will  contain  the  differences, 
in  model -surface  covariant  components,  between  a  projected  point  and  the 
unconstrained  L.S.  point  P.  It  can  be  shown  that  the  elements  in  the  guide 
vector  must  be  ^0.  otherwise  the  pertinent  branch  leads  to  a  non-L.S.  point. 

This  property  can  be  first  demonstrated  for  the  constrained  L.S.  point  P, 
which  can  be  thought  of  as  lying  in  a  u-m  dimensional  hyperplane.  Since  P  is 
the  result  of  an  appropriate  orthogonal  projection  (composed  of  m  ""closely 
nested"  projections)  of  P  onto  this  hyperplane,  a  u  dimensional  sphere  centered 
at  P  and  having  the  smallest  possible  radius  touches  the  hyperplane  at  P.  At 
this  stage,  we  Introduce  a  u  1  dimensional  hyperpJatu* ,  called  the  P  divider, 
which  is  tangent  to  this  sphere  at  P  and  generates  two  u-dimensional  half 
spaces.  The  P-divlder  itself  contains  the  u-m  dimensional  hyperplanc  of  the 
point  P. 

If  all  the  coordinate  axes  are  parallel  transported  from  the  origin  to  P, 
u-m  of  them  will  lie  in  the  above  u-m  <limensional  hyperplane  (always  containing 
the  axes  of  coordinates  not  subject  to  constraints),  and  the  remaining  m  must 
lie  in  the  half -space  excluding  the  point  P  But  this  means  that  the  cosine  of 
the  angles  formed  by  the  vector  PP  and  by  any  of  the  u  coordinate  axes  must  be 


either  0  (with  regard  to  the  u-m  axes  spanning  the  urn  dimensional  hyperplane  of 
P),  or  larger  than  zero  (with  regard  to  the  remaining  m  axes).  This,  in  turn, 
leads  to  u-m  zeros  and  to  m  values  larger  than  zero  in  the  corresponding  entries 
of  the  guide  vector  formed  for  the  constrained  L.S.  point  P. 

The  Intermediate  points  formed  by  the  projections  that  ultimately  generate 
P  can  be  thought  of  as  a  sequence  of  nested  L.S.  problems  with  certain 
constraints  removed.  Therefore,  the  entries  of  the  guide  vectcis  associated 
with  the  Intermediate  points  along  a  given  branch  must  be  >0,  otherwise  this 
branch  will  not  result  in  the  L.S.  point  P  and  should  be  rejected.  Since  a 
( u-dlmens ional )  sphere  touches  a  hyperplane  of  any  dimensions  in  one  point  only, 
there  can  be  only  one  L.S.  point  P.  in  conclusion,  if  any  branch  reaches  a 
point  for  which  all  components  subject  to  constraints  are  ^0,  and  for  which  all 
guide-vector  entries  ate  >0,  this  point  is  P.  The  unique  L.S.  solution  has 
thus  been  found  and  all  the  remaining  branches  should  be  discontinued. 


Description  of  the  Geometrical  Algorithm 

The  unrestricted  L.S.  solution  X  (I.e,,  a  column- vector  of  u  elements 
corresponding  to  the  point  P)  is  computed  from  (3a)  as 

X  -  il  .  (7) 


It  would  have  been  formally  more  appropriate  to  denote  the  unrestricted  L.S. 
solution  by  a  different  symbol,  eg.  X^^\  and  to  reserve  the  symbol  X  for  the 
constrained  L.S.  solution.  However,  the  geometrical  algorithm  will  be  described 
more  conveniently  with  the  unrest licted  solution  reprccsented  by  X,  and  with  the 
constrained  solution  attributed  superscripts  in  parentheses.  We  now  partition  X 
into  u-1  elements  and  into  the  remaining  one  element,  for  example  the  last,  and 
partition  the  column -vector  l'  in  the  same  way.  iml  I.-r  ’ y ,  ‘he  ">T*Ti.r  N  is 
partitioned  into  the  leading  submatrix  of  dimensions  (u  l)x(u-l)  and  into  the 
remaining  submatrices  of  dimensions,  clockwise,  (>i  1)<1,  1xi  (a  single  number), 
and  1x(u  ]).  These  partitions  arc  symbolized  by 
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(8a , b, c ) 


where  X  ,  L'  ,  and  N 

11  11  on 


arc  single  numbets. 


It  stu)uld  he  borne  in  wind  tiiat  tin?  suhs(;ript  "1"  does  not  ( orrespond  to 

one  element,  or  to  the  first  element,  but  to  .1  set  of  n-1  elements,  here  th«.‘ 

first  u  1  of  u  elements.  On  the  other  hand,  the  subscript  11"  pertains  to  one 

element,  here  the  last  of  n  elements.  Since  the  matrix  N  is  posi tive- def Inite , 

so  are  its  diagonal  submatrict-s  N..  and  N  ;  here,  in  particular.  N  >0.  The 

t 1  un  uu 

snbmatrices  and  are  transposes  of  each  olhei  A  computational  advantage 
stemming  from  an  arrangement  such  as  (8c)  is  that 
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(8d) 


the  reciprocal  value  of  a  /lumber.  Most  importantly,  partitions  of  this  kind 
correspond  to  the  geometrical  strategy  of  closely  nested  pr/)ject i ons .  The 
formulas  presented  in  matrix  notation  in  the  remainder  of  this  section  have  been 
developed  via  geometry,  although  some  of  them  can  be  derived,  in  a  more  tedious 
manner,  also  by  algebraic  means. 


If  the  unrestricted  I,.S.  solution  produces  X  <0.  tht*  constraint  X  -0  may 

u  u 

be  binding.  Henceforth  we  assume  the  constraints  (6b).  ot  (•‘3b)  where  X^  belongs 
to  the  set  X".  Thus,  a  reference  to  an  element  as  being  "negative"  implies 
"conflicting  negative".  Should  the  above  constraint  be  indeed  binding,  it  would 
Induce  a  projection  onto  the  appropriate  hyperplane  of  dimensions  u-1.  This 
possibility  should  be  scrutinized  in  an  appropriate  branch.  As  an  example,  In  a 
three-dimensional  model  hyperplane  with  oblique  Cartesian  axes  (x.y.z),  an 
unconstrained  result  z<0  (with  respect  to  as  the  origin)  would  induce  the 
projection  of  P  onto  the  plane  (x.y)  in  one  of  possible  branches. 


In  considering  as  the  last  parameter,  the  projection  of  P  onto  the 
hyperplane  spanned  by  the  first  u  1  coordinate  axes  can  be  shown  to  be 
equivalent  to  solving 

xf^’  (’  ,  X*^^  =  0  .  (9a.  b) 

1  111  u 

where,  in  analogy  to  {8a),  contains  u-1  elements  and  ^  is  a  single 

number.  The  variance-covariance  matrices  for  these  two  subsets  are  and 
0,  respectively.  However,  instead  of  Inverting  at  the  first  level,  and 
Inverting  progressively  smaller  matrices  at  further  levels,  one  can  take 
advantage  of  the  property  paralleling  (8d)  and  deduce 
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where  \  N**^,  and  are  snhmatiices  of  N 

De  partitioned  in  a  complete  analogy  to  (fir),  except 
replace  the  subscripts. 


the  lattei-  is  imagined  to 
that  the  superscripts  now 


We  next  attribute  to  X  a  supeiscript  in  parentheses,  indicating  tiie  level 
of  the  algorithm,  i.e.,  the  number  of  proj<!ctlons  effectuated  up  to  and 
including  the  current  step.  This  notation  has  already  been  used  in  (9a, b)  above 
for  the  first  level  of  projections.  In  general,  the  notation  xl*"^  identifies 
a  vector  of  u-m  elements  obtained  at  the  m  th  level:  the  remaining  m  elements  in 
the  complete  vector  X^™^  are  zero  The  elements  brought  to  zero  by  successive 
projections  are  not.  in  general,  the  last  elements  in  given  partitions  as  was 
assumed  in  (8a-d)  and  (9a, b)  for  convenience.  Instead,  they  are  the  negative 
elements  which  have  effectively  induced  such  projections.  However,  their 
location  in  X  or  its  subsequent  part  it  ions  does  not  alter  the  architecture  of 
formulas  such  as  (10).  We  only  have  to  make  sure  that  N*"'  corresponds  to  the 
negative  element  inducing  the  (next)  projection.  It  is  now  denoted  \  where 
the  generic  index  ’’ i "  symbolizes  the  negative  element  in  question  Clearly,  ^ 
may  correspond  to  the  last  element  as  well,  provided  the  latter  turns  out  to  be 
negative,  and  provided  we  are  treating  the  branch  where  this  eJemenl  actually 
Induces  a  projection.  This  reflects  the  possibility  that  there  may  exist  other 
negative  elements  in  X,  each  of  which  may  induce  a  legitimate  (non  repeti tious ) 
projection  and  thereby  create  a  separat(>  branch. 


The  submatrices  of  the  type  N  |  will  also  be  attributed  a  superscript 

in  parentheses,  likewise  indicating  the  level  of  projections.  Thus,  N  is 
( 1  ) 

to  be  potentially  partitioned  into  the  submatrices  denoted 


written  as  N 

N(i)n  ^(1)11 

.,(1) 


and  and  equation  (10)  is  rewritten  as 


.11 


,1  i 


i/n“  )  n‘  ' 


(11) 


The  lack  of  superscripts  in  parentheses  on  the  right  hand  side  of  (11)  is 
equivalent  to  a  superscript  ' (0)”.  This  means  that  no  modification  has  taken 


place  as  yet,  and,  accordingly,  that  N 

1 


1 1 


.1  i 


'  ,  and  are  submatrices  of 


the  original  matrix  N 


In  general .  we  have 


..(m  1  )  1  i 


l/^(m-l)iij  ^(m  Dil 


,^,(m  1)11 


(12) 


(m  1)11 


arc 


where  the  dimensions  of  N*"*'  are  (n-m)x(u  m). 
the  same,  while  the  dimensions  of  and 

and  lx(ii  m).  Predictably.  N*""  is  a  single 


The  dimensions  of  N 

„(m  1  )  i1  .  ,  , 

N  are  rc.sj»c(  1 1  vcl  y 

numbf^r  ri-(;ar;l  I  css  of  m. 


( u  m )  X  1 


Although  the  vector  IJ  has  served  in  obtaining  the  unrestricted  L.S. 
solution  in  (7),  it  will  not  be  used  in  any  capacity  for  any  other  task.  This 
stems  from  the  fact  that  the  geometrical  algorithm  avoids  new  matrix  inversions 
such  as  featured  in  (9a),  and  takes  advantage  instead  of  the  relations  such  as 
(12),  involving  relatively  very  few  scalar  multiplications.  In  lieu  of  {9a), 
the  first-level  solution  is  given  by  the  algorithm  as 

=  Xj  (1/N^^)  X.  ,  (13) 

and  in  lieu  of  inverting  one  can  obtain  the  variance  covariance  matrix 

associated  with  this  solution  as  in  (11).  In  general,  the  m-th  level  solution 
is  given  by  the  algorithm  as 

J,(»I  .  u  .  nil  ^(»  I)  1,^1 

( Bi ) 

with  the  corresponding  variance-covariance  matrix  N  '  presented  in  (12).  The 
remaining  m  element.^  in  X'  ^  are  zero,  and  the  pertinent  variance-covariance 
matrix  is  a  zero  matrix. 


We  notice  that  not  only  does  the  path  from  (11)  to  (12)  involve  relatively 
few  scalar  multiplications  as  has  already  been  mentioned,  but  these  operations 
take  place  in  successively  smaller  systems.  A  similar  statement  applies  also 
when  proceeding  from  xj^^  to  xj"**.  We  reiterate  that  the  generic  index  "i" 

In  all  of  the  above  formulas  identifies  a  negative  element  of  the  solution 
vector  at  a  given  level.  This  element  is  necessarily  different  from  level  to 
level.  In  fact,  once  a  negative  element  has  been  suppressed  by  a  projection,  it 
becomes  zero  and  remains  fixed  at  that  value. 


Without  the  aid  of  the  guide  vector,  we  would  attempt  to  reach  the 

coi. strained  L.S.  solution  subject  to  (6b),  for  example,  by  examining  various 

branches  while  avoiding  equivalent  (l.e.,  essentially  repetitious)  paths.  This 

would  lead  to  one  or  more  solutions  at  a  given  level,  e.g.  the  level  m,  the 

T 

elements  of  which  would  be  ail  non  negative.  The  values  of  AV  PV  due  to  each 

projection  along  an  individual  path  would  be  accumulated  and  added  to  the  value 
T 

of  V  PV  from  (2),  associated  with  the  unrestricted  adjustment.  The  decisions  as 
to  whether  a  given  outcome  is  a  L.S.  solution,  and  as  to  which  branches  can  be 


T 

discontinued  at  what  level,  would  then  be  made  based  on  the  accumulated  V  PV. 

T 

The  Increment  in  V  PV  due  to  the  m-th  projection  is  very  sl^mple  to  compute  (ii 
is,  in  fact,  already  partially  evaluated  in  equation  14): 


(15al 


T 

Based  on  (15a).  the  accumulated  V  PV  is  computed  as 


p  V  ^  ^ 


T  T  ( 1 

V  PV  *  AV  PV 


AV^PV^^’  .  .  ,  .  +  AV^PV*™* 


(15b) 


{  in  ) 

If  a  certain  branch  yielded  x'  ^>0.  all  the  branches  where  the  accumulated 
T 

value  of  V  PV  has  surpassed  the  current  value  (15b)  for  this  solution  would  be 

discarded.  The  remaining  branches  would  be  continued,  since  the  accumulated 
T 

V  PV  in  one  or  more  of  them  could  be  smaller  than  the  current  value  (even  though 
they  might  entail  more  than  m  projections).  If  the  latter  should  occur,  it 
would  indicate,  in  turn,  that  the  current  branch  has  resulted  in  a  non-L.S. 
solution.  In  spite  of  the  simplicity  of  (15a. b),  this  strategy  has  the  drawback 
that  it  discards  a  branch  "after  the  fact”,  i.e,,  after  an  additional  projection 
has  been  computed. 

In  this  context,  the  guide  vector  acts  as  an  invisible  arm  reaching  forward 
into  the  next  projection  and  detecting  a  possible  cause  for  rejection  (when  one 
or  more  of  the  vector's  elements  are  negative).  Suppose  that  the  current  m-th 
level  has  resulted  in  the  solution  xj'''^  obtained  upon  using  etc., 

according  to  (14).  Suppose  further  tliat  some  of  the  elements  in  are 

negative.  However,  Instead  of  a  complete  formulation  of  one  distinct  matrix 

associated  with  each  negative  element  in  view  of  computing  distinct  vectors 
jj(m+l)  giving  rise  to  distinct  branches,  we  compute  only  one  diagonal 
element  of  each  such  matrix  . 


For  the  sake  of  illustration,  we  consider  one  (possibly  the  only)  negative 
(  H)  ) 

element  of  X^  .  which  we  symbolize  by  "j",  i.e.,  we  consider  the  element 
Xj*^<0.  (The  negative  element  at  the  previous  level  was  Xj*  ^^.)  The 
current  negative  element  xj"'^  has  the  potential  to  generate  the  branch  "j"  at 
the  next  level,  numbered  m^l,  Ttiis  branch.  If  treated,  will  require  the 
knowledge  of  the  pertinent  matrix  in  order  to  compute  the  solution  vector 

Howevei-,  we  only  compute  the  j-th  diagonal  element  of  N*‘”^  wliich  is 
needed  as  a  building  f)lock  of  Mie  guide  vectoi'  at  this  stage.  In  denoting  it 
j^(m)jj  consulting  (12),  we  deduce  that 
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l)jj 


(16) 


(n‘®  ’  l^f  ‘  ‘  J 

The  remaining  elements  of  for  the  hi  ati<  h  "j"  will  bf-  Loroputeci  only  if  this 

branch  is  not  rejected  by  the  guide  vector. 

We  are  now  in  a  position  to  complete  the  formation  of  the  guide  vector.  As 
has  been  explained,  this  vector's  role  is  to  examine,  from  a  given  level,  a 
potential  projection  to  be  performed  at  the  next  level.  The  guide  vector  is  the 
most  easily  formed  at  the  level  1.  although  at  this  level  it  cannot  provide  any 
service  as  will  become  clear  presently.  In  assuming  that  the  element  "i"  in  the 
original  vector  X  is  negative,  Xj<0,  and  in  proceeding  along  this  branch,  one 
can  show  that  all  the  entries  of  the  guide  vector  at  the  first  level  are  zero 
except  for  the  entry  "i".  which  is  positive.  The  latter  is  obtained  as 

-  (1/N^‘)(-X.)  >  0  .  (17) 

Since  by  construction,  It  is  incapable  of  rejecting  its  bianch. 

We  confirm  the  above  assertion  for  a  three  dimensional  model  hyperplane 
with  oblique  Cartesian  axes  {x,y,z},  where  an  unrestricted  L.S.  point  P  has 
z<0.  The  projection  of  P  onto  the  plane  {x.y}  generates  the  point  P',  for  which 
the  first  two  covariant  components  of  the  guide  vector  are  zero.  Regardless  of 
whether  P'  falls  inside  the  desired  (sub-)region  x>0,  y>0,  the  corresponding 
P'-divider  is  Identical  with  the  plane  {x,y}.  while  the  z-axis  points  away  from 
the  half -space  of  P;  the  latter  point  is  "below"  the  plane  (x.y)  by  virtue  of 
z<0.  But  since  the  z-axis  and  the  guide  vector  are  pointing  into  the  same  half¬ 
space,  the  third  covariant  component  of  this  vector  is  positive.  Similarly.  P 
having  a  negative  i-th  coordinate  may  bo  projected  from  a  u -dimensional  model 
hyperplane  onto  a  u-1  dimensional  embedded  hyperplane,  generating  P'  as  well  as 
the  corresponding  P'  divider.  The  latter  coincides  with  the  embedded  hyperplane 
and  thus  contains  u-1  of  the  u  coordinate  axes,  while  the  remaining  axis  "1" 
points  away  from  the  half  space  of  P. 

We  next  suppose  that  the  element  "j"  of  xj ' *  is  negative,  j.e.,  Xj^^<0. 

Before  proceeding  to  the  full-scale  compulation  at  the  second  level  (branch  "j") 

necessitating  in  order  to  produce  X^^^.  we  compute  the  guide  vector  G 

for  such  a  solution  in  order  to  determine  whether  this  branch  should  not  be 

( 2 ) 

discarded  altogether.  The  vector  (>  will  have  two  entries  different  from  zero 
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("i”  and  "j"),  while  all  its  remaining  entries  will  be  zero.  The  entry  "j"  will 
be  positive,  computed  in  analogy  to  (17): 


’  1  >  0  .  (18) 

where  would  be  obtained  from  the  formula  (16)  with  m=l.  The  entry  "i"  is 

composed  of  two  parts,  namely  the  previous  entry  "i  "  plus  a  simple  correction; 

-  (1/N^^  .  (19) 

1  1  J 

( 2 ) 

If  G,  turns  out  to  be  negative  (due  to  the  correction),  the  branch  "i",  "j" 

^  ( 1 )  ( 2 ) 

is  rejected.  Otherwise  one  proceeds  to  complete  N  ,  leading  to  X 

( 2  ) 

If,  in  this  case,  an  element  "k”  of  the  solution  X  is  negative,  i.e.,  if 

(2)  ( 3 ) 

Xj^  <0,  we  generate  the  guide  vector  G  to  test  the  third  level.  This 

vector  will  have  three  nonzero  entries,  "i",  "j",  and  "k”.  The  entry  "k"  is 

positive,  computed  in  analogy  to  (17)  or  (18): 

G^^’  =  [l/N^^’‘^‘'][-xJ.^’ ]  >  0  ,  (20) 

K  K 

( 2  )kk 

where  N  would  be  readily  obtained  from  (16)  with  di=2  and  with  j  and  k 

replacing  1  and  j,  respectively.  The  entry  "J"  is  composed  of  the  previous 
entry  "j"  plus  a  simple  correction  (one  term): 

o'^'  .  of  ,  (211 

J  J  k 

Clearly,  the  elements  etc.,  are  available  from  the  previous 

level  (see  the  statements  following  equation  19).  Finally,  the  entry  "i"  is 
composed  of  the  previous  entry  "  i "  plUvS  a  correction  consisting  of  two  terms: 

g(3)  ^  ^(2)  c!^^]  .  (22) 

i  1  J  k 

It  Is  noteworthy  that  all  the  building  blocks  needed  in  (20)-(22)  are  available, 

( 2 )  kk 

whether  beforehand  or  sequentially,  with  the  exception  of  N  ,  which  requires 
a  very  simple  computation  as  indicated  above. 

The  hierarchy  in  building  the  guide  vector  at  higher  levels  is  now  quite 
apparent.  The  most  recent  entry  in  G^"*^  is  always  positive,  and  its  computation 
requires  an  element  obtainable  via  (16).  The  ne,xt  most  recent  (nonzero)  entry 
in  G^"^  is  the  corresponding  entry  in  G^"*  plus  one  correction  term.  The 
following  entry  in  G^™^  is  the  corresponding  entry  in  G^"  plus  two  correction 


terms,  etc.  The  correction  terms  in  the  last  nonzero  entry  of  G 


contain 


19 


elements  of  the  original  matrix  N  .  l.o.,  elements  belonging  to  the  level  "0". 
(  fi  ) 

[f  any  entry  In  G  is  negative  at  any  level,  £  =  2,3 . m .  the  pertinent 

branch  Is  rejected.  We  recall  that  if  any  branch  contains'  a  permutation  of  a 
sequence  of  projections  already  rejected,  such  a  branch  should  be  discarded 
without  testing. 


Numerical  Example  Solved  by  the  Geometrical  Algorithm 

The  simple  example  below  is  presented  at  the  level  of  normal  equations 
(3d),  where  the  inequality  constraints  have  the  form  (6b): 


N  X  =  U  . 
X  >  0  . 


(23a) 

(23b) 


The  number  of  parameters  is  four,  I.e.,  u-4.  The  matrix  N  and  the  vector  U  are 
given  as 


0.5  0  3.75 


N  =  0.5 

0 


-3.75  0 


0.675 

U  =  0.35 

0.2 

5.1 


(24a, b) 


The  matrix  N  is  positive-definite,  and  as  such  can  be  decomposed  into  the 
T 

product  T  T,  where  the  (real)  matrix  T  is  regular,  upper -triangular . 

The  unrestricted  L.S.  solution,  denoted  X  in  agreement  with  an  earlier 
convention,  follows  from  (23a): 

X  =  U  .  (2J 


In  particular,  we  have 


944 

432 

-40 

148 

-1 

-432 

444 

55 

56 

X  --  0 

-40 

-55 

75 

-18 

0 

148 

-56 

-18 

32 

•> 

-0 

itr  lx 

N  can 

bo  1 

obtained 

from  N  as  T 

0.35 


1  .T 


(26a,b) 


was  introduced  above;  T  ^  (real)  is  likewise  regular  and  upper  triangular. 


In  applying  tin-  geometrical  algorithm  of  the  previous  section,  we  replace 
the  compact  notation  x| ^ ^  by  a  more  explicit  notation 

when  we  need  to  Identify  specific  elements.  Thus,  in  the  branch  1  below,  k=2, 
8=3,  and  m=4.  This  kind  of  superscript  notation  will  be  used  at  any  level 

(including  the  level  0),  but  only  in  conjunction  with  various  vectors  "X"  and 

I  4 

their  subsets.  Since  X  -  -1.6  and  X  =0.5,  the  level  1  will  give  rise  to  two 
branches,  namely  branch  1  (named  after  the  element  xS  and  branch  4  (named  after 

4 

the  element  X  ).  Branch  1  will  correspond  to  the  possibly  binding  constraint 
X^=0,  or,  more  precisely,  X^^^^  0,  and  to  the  projection  of  the  unrestricted 
L.S  point  P  onto  the  three-dimensional  hyperplane  spanned  by  the  second,  third, 

and  fourth  coordinate  axes.  As  has  been  just  describt.-d,  the  vector  X^^^  is 

r„(l)2  „(1)3  „(1)4,T  •  rv2  „3  „4,T 

[X  ,  X  ,  X  ]  ,  while  the  corresponding  vector  X^  is  [X  ,  X  ,  X  ] 

=[0.8,  0.35,  0.5]^.  The  solution  x|'^  then  follows  from  (13),  where  i=l, 

i.e.,  where  X.  is  X  --1.6,  Similar  comments  with  self-evident  modifications 
1 

apply  also  for  the  branch  4,  where  i  4. 

If  the  branch  1  resulted  in  X^^^>0,  and  thus  X^^^>0  (due  to  X^^^^=0), 

the  unique  constrained  L.S.  solution  would  be  achieved.  In  this  case,  the  guide 

vector  composed  of  from  (17)  and  of  u-l  zeros  would  be  unnecessary 

(all  its  entries  would  be  >0).  However,  (17)  would  be  a  stepping  stone  in  the 
T 

computation  of  AV  PV  as  given  by  (15a),  since 
T  ( 1 )  ( 1 ) 

AV  PV'  '  =  G:  '  (  X.)  ,  (27) 

II 

where  X^  is  X^=  1.6  as  above,  with  similar  relationships  occurring  at  higher 
levels  as  well.  On  the  other  hand,  if  the  branch  1  does  not  lead  to  such  a 
quick  solution,  guide  vectors  for  the  next  level  should  be  formed  for  all 
possible  new  branches. 

Branch  1 .  In  accordance  with  the  above  description,  we  form 


-X(l)2- 

0.8 

-432 

0.067797 

X(l)3 

0.35 

(1/173) 

-40 

(173/944)(-l .6)  = 

0.282203 

x(l)4 

-0.5 

148 

-0.249153 

(1)4 

From  X  <0,  we  conclude  that  at  least  one  more  level  will  be  needed.  The  next 
level  (I.e.,  the  level  2)  will  give  rise  to  the  branch  1,4.  However,  before 
proceeding  any  further,  we  form  the  guide  vectors  0^'^  and 


the  latter  to 


test  whether  the  branch  1,4,  and  thus  the  entire  branch  1,  should  not  be 
rejected.  First,  (17)  yields 


=  ( 3  73/944  )x  1  f,  .  0.293220  . 

T  ( 1 )  12) 

from  which  it  follows  that  AV  PV  =0.469353.  In  view  of  G  ,  we  compute 
( 1  )44 

N  via  (16)  with  m=l.  upon  the  substitution  i  =  l  and  J  =  4 ; 

N^^*'*'*  =  (1/173)  (32  -  148^/944)  -  0.050847  . 

With  this  entry.  (18)  and  (19)  yield 

=  (1/0. 050847 )x0. 249153  =  4.9  , 

4 

=  0.293220  -  ( 1 73/944 ) ( 1 48/1 73 ) x 4 , 9  =  -0  175. 

Accordingly,  the  entire  branch  1  is  rejected. 

Branch  4.  This  branch  corresponds  to  the  possibly  binding  constraint 

( 1  )4 

X  =0.  In  analogy  to  the  preceding  procedure,  we  have 


-1.6 

148 

0.7125 

^(1)2 

= 

0.8 

-  (1/173) 

-56 

(173/32) (-0.5)  = 

-0.075 

0.35 

- 18 

0.06875 

The  only  possible  branch  at  the  level  2  Is  4,2.  We  notice  that  since  the  branch 
3,4  was  rejected  above,  so  would  be  the  branch  4,1;  but  the  latter  would  now 
correspond  to  "projecting  a  positive  component",  which  had  been  eliminated  from 
the  strategy.  In  view  of  the  guide  vectors  and  we  compute 

=  (173/32)x0.5  =  2.703125  , 

4 

from  which  it  follows  immediately  that  AV^PV* ^ ^ = 1 . 351562 .  Further,  with  m=l 
and  the  substitution  l--=4  and  j=2,  equation  (16)  yields 

=  (1/173)(444  56^^32)  =  2  . 

With  this  element,  It  follows  from  (18)  and  (19)  that 

=  (l/2)x0.075  0.0375  . 

g1^’  =  2.703125  -  (173/32)(-56/173)x0.0375  =  2.76875  . 

4 


Thus,  the  branch  4,2  at  the  level  2  will  not  be  rejected.  Furthermore,  upon 

i  2 1  T  <  2  ) 

considering  (15a),  the  above  result  for  '  yields  AV  PV'  =0.002813. 


Branch  4,2.  This  branch  corresponds  to  the  possibly  binding  constraint 

(2)2  ( 1 ) 

X  =0.  We  first  use  (11)  with  i=4  to  compute  N  .  which  will  serve  to 

12\ 

evaluate  x' 


,(1) 


(1/173){ 


944 

432 

-40 


432 

444 

-55 


-40 

-55 

75 


'148* 

56 

-18 

(1/32)[148  -56  -18]} 


N 


(1  ) 


1.5  - 1  0.25 

-1  2  -0.5 

0.25  -0.5  0.375 


This  result  Is  the  same  as  the  Inverse  of  the  (3x3)  leading  submatrix  of  N  in 
(24a).  However,  the  above  scheme  is  much  more  economical,  amounting  to  only  a 
few  scalar  multiplications  (advantage  should  be  taken  of  the  symmetry  of  N  \ 
etc.  )  . 

With  available,  the  solution  follows  from  (14)  as 


■x(2)l‘ 

0.7125 

-1 

0.675 

,(2)3 

L-  J 

0.06875j 

_ 

-0.5 

(l/2)(-0.075)  = 

0.05 

Since  these  values  are  positive,  the  constrained  L.S.  solution  has  been 
achieved,  namely 

X^^^  =  (0.675  0  0.05  OJ^  . 

T 

According  to  (15b),  the  accumulated  V  PV  is 

V^PV^^^  =  V^PV  +  1.35)562  ^  0.002813  =  V^PV  +  1.354375  . 


We  note  that  if  the  algorithm  started  with  the  branch  4.  branch  4,2  would  yield 
the  unique  constrained  L.S.  solution  ami  branch  1  would  be  skipped  altogether. 


IJescriiJt  i(Mi  of  l  hi^  Oua<lr«)(.i r  Pi  oKt  aMwi  ng  AlHoi  iHim 

The  I. .  S .  adjustment  of  a  linear'  pafamet  lii-  mndcl  with  linear-  inequality 
ronstraints  ean  be  treated  by  the  method  r>f  ijuaflrat  i c:  proj’.rammint;  presented  in 
[Bard,  1074],  Seetivin  6-2.  In  principle,  the  prrjblem  formulated  on  page  147 
therein  seeks  to  minimize  the  quadratic  function  Q  in  the  parameters  X,  namely 

Q(X)  =  N  X  -  2  U  c  .  (28a) 

subject  to  the  inequality  constraints 

C  X  »  0  .  (28b) 

The  dimensions  of  N.  C,  X,  and  I’  are  respectively  uxu,  s^u,  uxi,  and  ux l .  This 

T 

task  is  equivalent  to  minimizing  V  PV  as  in  (2),  subject  to  the  constraints 

T 

(28b).  Accordingly,  in  terms  of  our  notation,  r  in  (28a)  corresponds  to  L  PL. 

It  should  be  mentioned  that  the  notation  in  (Bard.  1974]  differs  from  that 
employed  herein  in  several  respects.  For-  e.xample,  v  corresponds  here  to  X,  R 
corresponds  to  N,  q  corresponds  to  tJ,  etc. 

As  in  (25),  we  denote  the  unrestricted  I..S.  solution  by  X  and  write 

X  =  n'^  1]  . 

Thr?  algorithm  presented  by  Bard  [1974]  is  transcr  ibed  below  in  otrr  notation. 

1)  Form  the  matrix  E  of  dimensions  sx{sf]): 

E  s  [W  z]  =  (CN“^  CX]  ;  (29) 

W  is  a  matrix  of  dimensions  sxs,  and  z  is  a  col umn -vector  of  s  elements. 

2)  In  conjunction  with  E,  form  the  vector  k  of  s  elements,  whosj'  enli ies  are 
initially  set  to  unity. 

3)  Find  the  "i"  for  which  z ,  k  .  -  a  =  m inimum ;  initially,  a  will  be  the  smallest 

i  1 

element  of  z.  If  a?0,  go  to  step  5. 

4)  if  a<0,  execute  the  Gauss  .Jordan  pivot  on  the  element  (i,i)  of  F, ,  and 
change  the  sign  of  k.  .  With  the  notation  t-  1/E.j,  this  pivot  trig  consists 
of  two  basic  steps,  where  symbolizes  the  replacement  of  elements: 

a)  For  all  the  elements  of  E  except  the  i  th  row  and  the  i  th  column; 

E  -  E  -  E  ,  E.  t  . 
pq  pq  pi  iq 
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bj  For  aJ]  the  romainirif;  «-Ienients,  start  inn  with  the  i-th  row  except 

the  eJement  (  I ,  i )  : 

K  .  »  E ,  t  ; 

tq  iq 

f.ont  irminH  with  the  i  th  column  except  the  eJement  (i.i): 

E  -  -E  .  t  : 

pi  pi 

and  ending  with  the  element  (i,i): 

E.  .  t  . 

1 1 

in  fad,  the  liuee  p<»r  t  c  of  the  current  iter.,  "b"  eould  be 

performed  In  any  order. 

The  new  matrix  E  and  the  new  vi^ctor  k  replace  their  previous  counterparts, 
and.  accordingly,  are  denoted  again  E  and  k.  Within  E,  the  symbols  W  and  z 
are  preserved  as  well.  After  the  above  replacements,  return  to  step  3. 

5)  Suppose  s'  of  the  s  elements  in  k  are  -1  (the  others  being  a).  Form  the 
matrix  C  of  dimensions  s'xu  from  C.  upon  grouping  the  rows  of  tlie  latter 
corresponding  to  the  elements  -1  in  k.  This  procedure  eliminates  the  rows 
of  C  corresponding  to  +1  in  k,  which  represent  nonbinding  constraints.  In 
the  same  way,  form  the  column- vecior  z’  of  s'  elements  from  z  (z  is  the 
last  column  in  the  most  recent  matrix  E).  The  constrained  L.S.  solution, 
denoted  X' ,  is  given  as 

X'  -  (U  C'^  ■/.')  .  (30) 


Numerical  Example  Solved  by  the  Quadratic-Programming  Algorithm 

The  example  used  in  this  illustration  is  the  one  solved  previously  by  the 
geometrical  algorithm.  Accordingly,  Cl,  s^u=4,  and,  from  (29), 

E  =  [W  z]  -  IN  '  X]  . 

The  values  of  N  ’  and  X  are  given  in  (26a, b).  Thus,  the  initial  matrix  R  and 


the  Initial  vector  k  are 


5.456647 

2.497110 

() .  2.3 1 2  1 4 

0,8.5.5491  ' 

1.6 

’  1 

E 

2,497110 

2.566474 

0.317919 

0. 32.3699  ! 

0  .  H 

k  - 

1 

-0.231214 

-0,317919 

0.433526 

0.104046  I 

1)  .3,5 

1 

0,85.5491 

-0,323699 

0. 104046 

0.184971  ! 

-0.5 

1 

From  this  setup,  the  const ruiiuHl  L.S  soliiMnii  will  tx-  re'uchcd  in  four 
iterations,  in  the  sense  that  the  computations  will  pass  four  times  through  step 
3  of  the  quadratic  programming  algorithm. 

Iteration  1 .  The  Gauss  Jordan  pivot  will  be  executed  on  the  above  element 
(1  1).  for  which  t  =  1 /5 . 456647  - 0 . 183263 .  The  matrix  E  and  the  vector  k  become 


0.183263 

0.457627 

0.042373 

0  156780  ; 

0.293220 

-1 

E 

0.457627 

1,423729 

0.423729 

0.067797  , 

0.067797 

k  - 

1 

0.042373 

0.423729 

0.423729 

0.067797  1 

0  282203 

1 

-0.156780 

L 

0.067797 

0  067797 

0.050847  1 

-0  249153 

1 

This  iteration  corresponds  to  the  branch  1  of  the  geometrical  solution. 

Iteration  2.  The  pivoting  will  take  place  on  the  above  element  (4,4),  for 
which  t = 1 /O. 050847= 19.666667. 


0.666667 

-0.666667 

0. 166667 

3,083333  i 

J 

0.475' 

‘-1  ’ 

E  = 

0.666667 

1 .333333 

0.333333 

■1  333333  1 

0  4 

.  k  = 

1 

-0.166667 

-0.333333 

0.333333 

1 . 333333  1 

0.05 

1 

-3.083333 

1 .333333 

1.333333 

19.666667  1 

4.9 

-1 

This  iteration  would  correspond  to  the  branch  1.4  of  the  geometrnal  s('’utlon. 
if  this  branch  had  been  executed.  We  note  that  the  geometriiral  algorithm  would 
have  treated  a  smaller  system  at  this  stage,  namtrly  a  3x3  system.  By  contrast, 
the  quadratic-programming  algorithm  treats  essentially  a  4x4  system  at  every 
iteration.  We  have  seen  that  the  guide  vector  allows  the  geometrical  algorithm 
to  avoid  the  (wasteful)  execution  of  the  branch  1,4,  and  directs  it  instead  to 
the  remaining  branch  4. 

Iteration  3.  The  pivoting  will  take  plat-e  on  the  abovt-  element  (1,1),  for 
which  t  =  1 /0 . 666667 -- 1 . 5  .  This  is  the  second  time  the  element  (1.1)  is  subject 
to  the  Gauss  Jordan  pivot. 
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1  .5 

-  1 

0 . 25 

4.625  ' 

{ 

0.7125 

1 

E 

1 

2 

0.5 

1.75  1 

0.075 

k  = 

1 

0.25 

0.5 

0 . 375 

0 . 5625  1 

0.06875 

1 

4.625 

-1  .75 

0 . 5625 

5.40625  1 

2.703125 

-1 

This  Iteration  corresponds  to  the  t>ranch  4  of  the  ^jeonietr icol  solution. 

Iteration  4.  The  pivot i tig  in  this  last  iteration  will  take  place  on  the 
above  element  (2,2),  for  which  t  05. 


1 

0.5 

0 

3.75  1 

0.675 

= 

-  0 . 5 

0.5 

0.25 

0.875  1 

0.0375 

0 

0.25 

0.25 

1 

I  , 

0 . 05 

3.75 

0.875 

- 1 

6.9375  i 

-2 . 76875 

This  iteration  corresponds  to  the  branch  4,2  tif  the  geometrical  solution,  which 
yielded  the  constrained  L.S.  solution.  From  the  above  E  and  k,  the  quadratic- 
programming  algorithm  stipulat(*s  that 


'o  1 

0 

0.0375 

- 

,  Z '  ^ 

_0  0 

0 

Ij 

2.76875_ 

In  using  U  from  (24b),  we  obtain 

fl  C'*"  ■/.'  -  [0.675  ((.3875  0.2  2.33125]^  . 

Finally,  (30)  yields 

X'  -  n'*  ((!  -  C'^  z’)  fO.675  0  0.05  0)^  , 

whose  first  and  third  elements  have  already  appeared  at  their  respective  places 

in  the  last  column  of  E  above.  This  solution  X'  agrees  witli  the  final  result 
(2) 

X  obtained  in  the  branch  4,2  of  the  geometrical  solution. 
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